# State Minimum Wage and Mental Health Among Children and Adolescents
# Analyses using the Youth Risk Behavior Surveillance System
# N.M. Kavanagh, M. McConnell, N. Slopen
# September 20, 2024

# Please direct questions about this script file to nolankavanagh@fas.harvard.edu.

# Clear R environment
rm(list = ls())

# Load packages
library(here)         # Working directory
library(haven)        # Dataset tools
library(tidyverse)    # Analysis tools
library(dplyr)        # Analysis tools
library(psych)        # Analysis tools
library(lfe)          # Modeling tools
library(survey)       # Survey tools
library(ggplot2)      # Graphing tools
library(usmap)        # Graphing tools
library(cowplot)      # Graphing tools
library(scales)       # Graphing tools
library(arsenal)      # Table tools
library(modelsummary) # Table tools
library(kableExtra)   # Table tools
library(cdlTools)     # FIPS tools

##############################################################################
# Dataset preparation
##############################################################################

# Read YRBS datasets into R
yrbs_st_ad <- read_sav("YRBS data/SADC_2021_State_a_d.sav")
yrbs_st_eh <- read_sav("YRBS data/SADC_2021_State_e_h.sav")
yrbs_st_il <- read_sav("YRBS data/SADC_2021_State_I_L.sav")
yrbs_st_mp <- read_sav("YRBS data/SADC_2021_State_M_P.sav")
yrbs_st_qt <- read_sav("YRBS data/SADC_2021_State_Q_T.sav")
yrbs_st_uz <- read_sav("YRBS data/SADC_2021_State_U_Z.sav")

# Recode Arizona abbreviation
# Only relevant for A-to-D dataset
yrbs_st_ad <- yrbs_st_ad %>% mutate(
  sitecode_2 = case_when(
    sitecode == "AZB" ~ "AZ",
    TRUE ~ sitecode
  ))

# Convert state abbreviations to FIPS
# Necessary for a smooth merge
yrbs_st_ad$fipsst <- cdlTools::fips(yrbs_st_ad$sitecode_2, to="FIPS")
yrbs_st_eh$fipsst <- cdlTools::fips(yrbs_st_eh$sitecode,   to="FIPS")
yrbs_st_il$fipsst <- cdlTools::fips(yrbs_st_il$sitecode,   to="FIPS")
yrbs_st_mp$fipsst <- cdlTools::fips(yrbs_st_mp$sitecode,   to="FIPS")
yrbs_st_qt$fipsst <- cdlTools::fips(yrbs_st_qt$sitecode,   to="FIPS")
yrbs_st_uz$fipsst <- cdlTools::fips(yrbs_st_uz$sitecode,   to="FIPS")

# Bind all states together
yrbs <- bind_rows(yrbs_st_ad, yrbs_st_eh, yrbs_st_il, yrbs_st_mp, yrbs_st_qt, yrbs_st_uz)

##############################################################################
# Minimum wage dataset preparation
##############################################################################

# Read minimum wage dataset into R
min_wage_df <- read.csv("Supplemental data/State minimum wages, raw.csv")

# Note: The cleaning steps below produce virtually identical nominal wages for 1980-2019 as this source:
# https://raw.githubusercontent.com/Lislejoem/Minimum-Wage-by-State-1968-to-2020/master/Minimum%20Wage%20Data.csv
# Lisle did not have official data for 2020 onward, so we must clean the data ourselves.

# Copy forward missing years
# Note: These are only relevant for supplemental analyses
# We have complete information for the main study years (2000 onward)
min_wage_df$X1982 <- min_wage_df$X1981
min_wage_df$X1983 <- min_wage_df$X1981
min_wage_df$X1984 <- min_wage_df$X1981
min_wage_df$X1985 <- min_wage_df$X1981
min_wage_df$X1986 <- min_wage_df$X1981
min_wage_df$X1987 <- min_wage_df$X1981
min_wage_df$X1989 <- min_wage_df$X1988
min_wage_df$X1990 <- min_wage_df$X1988
min_wage_df$X1993 <- min_wage_df$X1992
min_wage_df$X1995 <- min_wage_df$X1994
min_wage_df$X1999 <- min_wage_df$X1998

# Convert dataset to long format
min_wage_long <- min_wage_df %>%
  select(order(colnames(min_wage_df))) %>%
  select(State.or.other.jurisdiction, X1980:X2023) %>%
  pivot_longer(cols = X1980:X2023, names_to = "year", values_to = "bls_min_wage")

# Convert state names to FIPS
# Necessary for a smooth merge
min_wage_long$fipsst <- cdlTools::fips(min_wage_long$State.or.other.jurisdiction, to="FIPS")

# Treat year as numeric
# Remove X before each value
min_wage_long$year <- as.numeric(gsub("X", "", min_wage_long$year))

# Generate new column for cleaned minimum wages
min_wage_long$state_min_wage <- min_wage_long$bls_min_wage

# Correct typo (double period) in select rows
min_wage_long$state_min_wage <- gsub("\\.\\.", "\\.", min_wage_long$state_min_wage)

# For ranges, select both values
min_wage_long <- min_wage_long %>% mutate(
  
  # Get values before hyphen or other symbols
  state_min_wage_1 = case_when(
    grepl("[-–/&]", state_min_wage) ~ gsub("[-–/&].*", "", min_wage_long$state_min_wage),
    TRUE ~ gsub("[-–/&].*", "", min_wage_long$state_min_wage)),
  
  # Get values after hyphen or other symbols
  state_min_wage_2 = case_when(
    grepl("[-–/&]", state_min_wage) ~ gsub(".*[-–/&]", "", min_wage_long$state_min_wage),
    TRUE ~ gsub(".*[-–/&]", "", min_wage_long$state_min_wage))
)

# Remove any remaining symbols or letters
# In practice, discard text starting with symbol
min_wage_long$state_min_wage_1 <- gsub("[-–&\\[\\(].*$", "", min_wage_long$state_min_wage_1)
min_wage_long$state_min_wage_2 <- gsub("[-–&\\[\\(].*$", "", min_wage_long$state_min_wage_2)

# Remove leading dollar signs and treat as numeric
min_wage_long$state_min_wage_1 <- as.numeric(gsub("\\$", "", min_wage_long$state_min_wage_1))
min_wage_long$state_min_wage_2 <- as.numeric(gsub("\\$", "", min_wage_long$state_min_wage_2))

# Get upper and lower ends of ranges
min_wage_long$state_min_wage_u <- pmax(min_wage_long$state_min_wage_1, min_wage_long$state_min_wage_2, na.rm = TRUE)
min_wage_long$state_min_wage_l <- pmin(min_wage_long$state_min_wage_1, min_wage_long$state_min_wage_2, na.rm = TRUE)

# Get federal minimum wages
min_wage_fed <- subset(min_wage_long, State.or.other.jurisdiction == "Federal (FLSA)")
min_wage_fed$federal_min_wage <- min_wage_fed$state_min_wage_1
min_wage_fed <- min_wage_fed %>% select(year, federal_min_wage)

# Merge state and federal minimum wages
min_wage_long <- left_join(min_wage_long, min_wage_fed, by="year")

# Get effective minimum wages
# Higher of state or federal minimums
min_wage_long$effective_min_wage_u <- pmax(min_wage_long$state_min_wage_u, min_wage_long$federal_min_wage, na.rm = TRUE)
min_wage_long$effective_min_wage_l <- pmax(min_wage_long$state_min_wage_l, min_wage_long$federal_min_wage, na.rm = TRUE)

# Read CPI dataset into R
# Source: https://www.usinflationcalculator.com/inflation/consumer-price-index-and-annual-percent-changes-from-1913-to-2008/
cpi_df <- read.csv("Supplemental data/CPI by month, raw.csv")

# Select January CPI estimates
# Note: Wage data are effective as of January 1 in each year
cpi_df$year        <- cpi_df$Year
cpi_df$cpi_january <- cpi_df$Jan
cpi_df <- cpi_df %>% select(cpi_january, year)

# Merge CPI estimates
min_wage_long <- left_join(min_wage_long, cpi_df, by="year")

# Compute inflation-adjusted wages
# Use CPI in January 2020: 257.971
min_wage_long$inflation_min_wage_u <- 257.971 / min_wage_long$cpi_january * min_wage_long$effective_min_wage_u
min_wage_long$inflation_min_wage_l <- 257.971 / min_wage_long$cpi_january * min_wage_long$effective_min_wage_l

# Reorder dataset
min_wage_long <- min_wage_long %>% arrange(fipsst, year)

# Generate lagged minimum wages
min_wage_long <- min_wage_long %>%
  group_by(fipsst) %>%
  mutate(
    
    # Not inflation adjusted
    lag_by_1  = lag(effective_min_wage_l, n=1, default=NA),
    lag_by_2  = lag(effective_min_wage_l, n=2, default=NA),
    lag_by_3  = lag(effective_min_wage_l, n=3,  default=NA),
    lag_by_4  = lag(effective_min_wage_l, n=4,  default=NA),
    lag_by_5  = lag(effective_min_wage_l, n=5,  default=NA),
    lag_by_6  = lag(effective_min_wage_l, n=6,  default=NA),
    lag_by_7  = lag(effective_min_wage_l, n=7,  default=NA),
    lag_by_8  = lag(effective_min_wage_l, n=8,  default=NA),
    lag_by_9  = lag(effective_min_wage_l, n=9,  default=NA),
    lag_by_10 = lag(effective_min_wage_l, n=10, default=NA),
    lag_by_11 = lag(effective_min_wage_l, n=11, default=NA),
    lag_by_12 = lag(effective_min_wage_l, n=12, default=NA),
    lag_by_13 = lag(effective_min_wage_l, n=13, default=NA),
    lag_by_14 = lag(effective_min_wage_l, n=14, default=NA),
    lag_by_15 = lag(effective_min_wage_l, n=15, default=NA),
    lag_by_16 = lag(effective_min_wage_l, n=16, default=NA),
    lag_by_17 = lag(effective_min_wage_l, n=17, default=NA),
    lag_by_18 = lag(effective_min_wage_l, n=18, default=NA),
    
    # Inflation adjusted
    lag_by_1_inf  = lag(inflation_min_wage_l, n=1,  default=NA),
    lag_by_2_inf  = lag(inflation_min_wage_l, n=2,  default=NA),
    lag_by_3_inf  = lag(inflation_min_wage_l, n=3,  default=NA),
    lag_by_4_inf  = lag(inflation_min_wage_l, n=4,  default=NA),
    lag_by_5_inf  = lag(inflation_min_wage_l, n=5,  default=NA),
    lag_by_6_inf  = lag(inflation_min_wage_l, n=6,  default=NA),
    lag_by_7_inf  = lag(inflation_min_wage_l, n=7,  default=NA),
    lag_by_8_inf  = lag(inflation_min_wage_l, n=8,  default=NA),
    lag_by_9_inf  = lag(inflation_min_wage_l, n=9,  default=NA),
    lag_by_10_inf = lag(inflation_min_wage_l, n=10, default=NA),
    lag_by_11_inf = lag(inflation_min_wage_l, n=11, default=NA),
    lag_by_12_inf = lag(inflation_min_wage_l, n=12, default=NA),
    lag_by_13_inf = lag(inflation_min_wage_l, n=13, default=NA),
    lag_by_14_inf = lag(inflation_min_wage_l, n=14, default=NA),
    lag_by_15_inf = lag(inflation_min_wage_l, n=15, default=NA),
    lag_by_16_inf = lag(inflation_min_wage_l, n=16, default=NA),
    lag_by_17_inf = lag(inflation_min_wage_l, n=17, default=NA),
    lag_by_18_inf = lag(inflation_min_wage_l, n=18, default=NA)
  )

# Generate indicator for above federal minimum
min_wage_long <- min_wage_long %>% mutate(
  above_fed = case_when(
    state_min_wage_l > federal_min_wage ~ 1,
    TRUE ~ 0
  ))

# Merge minimum wage and NSCH data
yrbs_all <- left_join(yrbs, min_wage_long, by=c("fipsst", "year"))

# Clear old datasets from R
rm(yrbs_st_ad, yrbs_st_eh, yrbs_st_il, yrbs_st_mp, yrbs_st_qt, yrbs_st_uz)

##############################################################################
# Supplemental Medicaid dataset preparation
##############################################################################

# Read Medicaid datasets into R
medicaid_1_5_df  <- read.csv("Supplemental data/Medicaid eligibility, ages 1-5, state-year, raw.csv")
medicaid_6_18_df <- read.csv("Supplemental data/Medicaid eligibility, ages 6-18, state-year, raw.csv")

# Add columns for 2001 and 2007
# Duplicate data from 2000 and 2006
medicaid_1_5_df$X2001  <- medicaid_1_5_df$X2000
medicaid_1_5_df$X2007  <- medicaid_1_5_df$X2006
medicaid_6_18_df$X2001 <- medicaid_6_18_df$X2000
medicaid_6_18_df$X2007 <- medicaid_6_18_df$X2006

# Restructure data to long format
library(reshape2)
medicaid_1_5_df  <- melt(medicaid_1_5_df, id.vars = c("X"),
                         variable.name = "year", value.name = "elig_1_5")
medicaid_6_18_df <- melt(medicaid_6_18_df, id.vars = c("X"),
                         variable.name = "year", value.name = "elig_6_18")

# Clean year variable
# Drop leading "X" in strings
medicaid_1_5_df$year  <- substring(medicaid_1_5_df$year,  2)
medicaid_6_18_df$year <- substring(medicaid_6_18_df$year, 2)

# Rename columns
colnames(medicaid_1_5_df)  <- c("state", "year", "elig_1_5")
colnames(medicaid_6_18_df) <- c("state", "year", "elig_6_18")

# Convert names to FIPS
# Necessary for smooth merge
medicaid_1_5_df$fipsst  <- cdlTools::fips(medicaid_1_5_df$state,  to="FIPS")
medicaid_6_18_df$fipsst <- cdlTools::fips(medicaid_6_18_df$state, to="FIPS")

# Treat year as numeric
medicaid_1_5_df$year  <- as.numeric(medicaid_1_5_df$year)
medicaid_6_18_df$year <- as.numeric(medicaid_6_18_df$year)

# Merge Medicaid and YRBS data
yrbs_all <- left_join(yrbs_all, medicaid_1_5_df,  by=c("fipsst", "year"))
yrbs_all <- left_join(yrbs_all, medicaid_6_18_df, by=c("fipsst", "year"))

##############################################################################
# Supplemental EITC dataset preparation
##############################################################################

# Read EITC dataset into R
eitc <- read.csv("Supplemental data/EITC policies, state-year, clean.csv")

# Merge EITC and YRBS data
yrbs_all <- left_join(yrbs_all, eitc, by=c("fipsst", "year"))

##############################################################################
# Supplemental TANF dataset preparation
##############################################################################

# Read TANF dataset into R
tanf <- read.csv("Supplemental data/TANF benefits family 3, state-year, clean.csv")

# Merge TANF and YRBS data
yrbs_all <- left_join(yrbs_all, tanf, by=c("fipsst", "year"))

##############################################################################
# Variable preparation
##############################################################################

# Treat fixed effects as factors
yrbs_all$year   <- as.factor(yrbs_all$year)
yrbs_all$fipsst <- as.factor(yrbs_all$fipsst)

# Adolescent's age
yrbs_all <- yrbs_all %>% mutate(
  age_2 = case_when(
    age == 1 ~ "12 years old or younger",
    age == 2 ~ "13 years old",
    age == 3 ~ "14 years old",
    age == 4 ~ "15 years old",
    age == 5 ~ "16 years old",
    age == 6 ~ "17 years old",
    age == 7 ~ "18 years old or older"
  ))

# Birth year
yrbs_all <- yrbs_all %>% mutate(
  birth_year = case_when(
    age == 1 ~ as.numeric(as.character(year)) - 12,
    age == 2 ~ as.numeric(as.character(year)) - 13,
    age == 3 ~ as.numeric(as.character(year)) - 14,
    age == 4 ~ as.numeric(as.character(year)) - 15,
    age == 5 ~ as.numeric(as.character(year)) - 16,
    age == 6 ~ as.numeric(as.character(year)) - 17,
    age == 7 ~ as.numeric(as.character(year)) - 18,
  ))

# Adolescent's sex
yrbs_all <- yrbs_all %>% mutate(
  sex_2 = case_when(
    sex == 1 ~ "Female",
    sex == 2 ~ "Male",
    TRUE ~ "Not provided"
  ))
yrbs_all$sex_2 <- factor(yrbs_all$sex_2, levels = c("Male", "Female", "Not provided"))

# Adolescent's race/ethnicity
yrbs_all <- yrbs_all %>% mutate(
  race_eth_2 = case_when(
    race7 == 1        ~ "American Indian or Alaska Native",
    race7 %in% c(2,5) ~ "Asian, Native Hawaiian, or Pacific Islander",
    race7 == 3        ~ "Black or African American",
    race7 == 4        ~ "Hispanic/Latino",
    race7 == 6        ~ "White",
    race7 == 7        ~ "Multiple races, non-Hispanic",
    TRUE ~ "Not provided"
  ))
yrbs_all$race_eth_2 <- factor(yrbs_all$race_eth_2, levels = c("American Indian or Alaska Native", "Asian, Native Hawaiian, or Pacific Islander", "Black or African American", "Hispanic/Latino", "White", "Multiple races, non-Hispanic", "Not provided"))

# Dichotomoize race/ethnicity
# Black/Latino vs. other for interaction models
yrbs_all <- yrbs_all %>% mutate(
  race_eth_cat = case_when(
    race7 %in% c(3:4)     ~ 0, # Black or Hispanic/Latino
    race7 %in% c(1:2,5:7) ~ 1  # Other races
  ))

# Adolescent's educational grade
yrbs_all <- yrbs_all %>% mutate(
  grade_2 = case_when(
    grade == 1 ~ "9th grade",
    grade == 2 ~ "10th grade",
    grade == 3 ~ "11th grade",
    grade == 4 ~ "12th grade",
    TRUE ~ "Not provided"
  ))
yrbs_all$grade_2 <- factor(yrbs_all$grade_2, levels = c("9th grade", "10th grade", "11th grade", "12th grade", "Not provided"))

# Sad or hopeless for 2+ weeks
yrbs_all <- yrbs_all %>% mutate(
  sad_hopeless = case_when(
    q25 == "1" ~ 1,
    q25 == "2" ~ 0
  ))

# Seriously considered suicide
yrbs_all <- yrbs_all %>% mutate(
  cons_suicide = case_when(
    q26 == "1" ~ 1,
    q26 == "2" ~ 0
  ))

# Made suicide attempt
yrbs_all <- yrbs_all %>% mutate(
  suicide_att = case_when(
    q28 %in% c("2":"5") ~ 1,
    q28 == "1"          ~ 0
  ))

# Alcohol use
yrbs_all <- yrbs_all %>% mutate(
  alcohol = case_when(
    q41 %in% c("2":"7") ~ 1,
    q41 == "1"        ~ 0
  ))

# Marijuana use
yrbs_all <- yrbs_all %>% mutate(
  marijuana = case_when(
    q47 %in% c("2":"6") ~ 1,
    q47 == "1"          ~ 0
  ))

# Physical fight
yrbs_all <- yrbs_all %>% mutate(
  fight = case_when(
    q16 %in% c("2":"8") ~ 1,
    q16 == "1"          ~ 0
  ))

# Lifetime minimum wage
# With and without inflation adjustments
yrbs_all <- yrbs_all %>% mutate(
  wage_life_nom = case_when(
    
    # 12 years or younger
    age == 1 ~ (effective_min_wage_l + lag_by_1 + lag_by_2 + lag_by_3 +
                  lag_by_4 + lag_by_5 + lag_by_6  + lag_by_7 +
                  lag_by_8 + lag_by_9 + lag_by_10 + lag_by_11 +
                  lag_by_12)/13,
    
    # 13 years
    age == 2 ~ (effective_min_wage_l + lag_by_1 + lag_by_2 + lag_by_3 +
                  lag_by_4  + lag_by_5 + lag_by_6  + lag_by_7 +
                  lag_by_8  + lag_by_9 + lag_by_10 + lag_by_11 +
                  lag_by_12 + lag_by_13)/14,
    
    # 14 years
    age == 3 ~ (effective_min_wage_l + lag_by_1 + lag_by_2 + lag_by_3 + 
                  lag_by_4  + lag_by_5  + lag_by_6  + lag_by_7 +
                  lag_by_8  + lag_by_9  + lag_by_10 + lag_by_11 +
                  lag_by_12 + lag_by_13 + lag_by_14)/15,
    
    # 15 years
    age == 4 ~ (effective_min_wage_l + lag_by_1 + lag_by_2 + lag_by_3 +
                  lag_by_4  + lag_by_5  + lag_by_6  + lag_by_7 +
                  lag_by_8  + lag_by_9  + lag_by_10 + lag_by_11 +
                  lag_by_12 + lag_by_13 + lag_by_14 + lag_by_15)/16,
    
    # 16 years
    age == 5 ~ (effective_min_wage_l + lag_by_1 + lag_by_2 + lag_by_3 +
                  lag_by_4  + lag_by_5 +  lag_by_6 +  lag_by_7 +
                  lag_by_8  + lag_by_9  + lag_by_10 + lag_by_11 +
                  lag_by_12 + lag_by_13 + lag_by_14 + lag_by_15 +
                  lag_by_16)/17,
    
    # 17 years
    age == 6 ~ (effective_min_wage_l + lag_by_1 + lag_by_2 + lag_by_3 +
                  lag_by_4  + lag_by_5  + lag_by_6  + lag_by_7 +
                  lag_by_8  + lag_by_9  + lag_by_10 + lag_by_11 +
                  lag_by_12 + lag_by_13 + lag_by_14 + lag_by_15 +
                  lag_by_16 + lag_by_17)/18,
    
    # 18 years or older
    age == 7 ~ (effective_min_wage_l + lag_by_1 + lag_by_2 + lag_by_3 +
                  lag_by_4  + lag_by_5  + lag_by_6  + lag_by_7 +
                  lag_by_8  + lag_by_9  + lag_by_10 + lag_by_11 +
                  lag_by_12 + lag_by_13 + lag_by_14 + lag_by_15 +
                  lag_by_16 + lag_by_17 + lag_by_18)/19
  ))
yrbs_all <- yrbs_all %>% mutate(
  wage_life_inf = case_when(
    
    # 12 years or younger
    age == 1 ~ (inflation_min_wage_l + lag_by_1_inf + lag_by_2_inf + lag_by_3_inf +
                  lag_by_4_inf + lag_by_5_inf + lag_by_6_inf  + lag_by_7_inf +
                  lag_by_8_inf + lag_by_9_inf + lag_by_10_inf + lag_by_11_inf +
                  lag_by_12_inf)/13,
    
    # 13 years
    age == 2 ~ (inflation_min_wage_l + lag_by_1_inf + lag_by_2_inf + lag_by_3_inf +
                  lag_by_4_inf  + lag_by_5_inf + lag_by_6_inf  + lag_by_7_inf +
                  lag_by_8_inf  + lag_by_9_inf + lag_by_10_inf + lag_by_11_inf +
                  lag_by_12_inf + lag_by_13_inf)/14,
    
    # 14 years
    age == 3 ~ (inflation_min_wage_l + lag_by_1_inf + lag_by_2_inf + lag_by_3_inf + 
                  lag_by_4_inf  + lag_by_5_inf  + lag_by_6_inf  + lag_by_7_inf +
                  lag_by_8_inf  + lag_by_9_inf  + lag_by_10_inf + lag_by_11_inf +
                  lag_by_12_inf + lag_by_13_inf + lag_by_14_inf)/15,
    
    # 15 years
    age == 4 ~ (inflation_min_wage_l + lag_by_1_inf + lag_by_2_inf + lag_by_3_inf +
                  lag_by_4_inf  + lag_by_5_inf  + lag_by_6_inf  + lag_by_7_inf +
                  lag_by_8_inf  + lag_by_9_inf  + lag_by_10_inf + lag_by_11_inf +
                  lag_by_12_inf + lag_by_13_inf + lag_by_14_inf + lag_by_15_inf)/16,
    
    # 16 years
    age == 5 ~ (inflation_min_wage_l + lag_by_1_inf + lag_by_2_inf + lag_by_3_inf +
                  lag_by_4_inf  + lag_by_5_inf +  lag_by_6_inf +  lag_by_7_inf +
                  lag_by_8_inf  + lag_by_9_inf  + lag_by_10_inf + lag_by_11_inf +
                  lag_by_12_inf + lag_by_13_inf + lag_by_14_inf + lag_by_15_inf +
                  lag_by_16_inf)/17,
    
    # 17 years
    age == 6 ~ (inflation_min_wage_l + lag_by_1_inf + lag_by_2_inf + lag_by_3_inf +
                  lag_by_4_inf  + lag_by_5_inf  + lag_by_6_inf  + lag_by_7_inf +
                  lag_by_8_inf  + lag_by_9_inf  + lag_by_10_inf + lag_by_11_inf +
                  lag_by_12_inf + lag_by_13_inf + lag_by_14_inf + lag_by_15_inf +
                  lag_by_16_inf + lag_by_17_inf)/18,
    
    # 18 years or older
    age == 7 ~ (inflation_min_wage_l + lag_by_1_inf + lag_by_2_inf + lag_by_3_inf +
                  lag_by_4_inf  + lag_by_5_inf  + lag_by_6_inf  + lag_by_7_inf +
                  lag_by_8_inf  + lag_by_9_inf  + lag_by_10_inf + lag_by_11_inf +
                  lag_by_12_inf + lag_by_13_inf + lag_by_14_inf + lag_by_15_inf +
                  lag_by_16_inf + lag_by_17_inf + lag_by_18_inf)/19,
  ))

# State has EITC program
yrbs_all <- yrbs_all %>% mutate(
  has_eitc = case_when(
    federal_pct > 0 ~ 1,
    TRUE ~ 0
  ))

# Generate nested sampling clusters
yrbs_all$cluster <- paste0(yrbs_all$fipsst, "-", yrbs_all$PSU, "-", yrbs_all$stratum)

##############################################################################
# Table: Demographic characteristics
##############################################################################

# Code inclusion/exclusion criteria
# That is, must have age and at least once outcome
yrbs_all <- yrbs_all %>% mutate(
  included_in_study = case_when(
    !is.na(age_2) & (!is.na(sad_hopeless) | !is.na(cons_suicide) | !is.na(suicide_att) |
      !is.na(fight) | !is.na(alcohol) | !is.na(marijuana)) ~ 1
  ))

# Assess missingness due to inclusion/exclusion criteria
# Note: Must subset to study years, i.e. 2001 to 2021
prop.table(table(subset(yrbs_all, year %in% c(2001:2021))$included_in_study, useNA="always"))

# Complete cases for models
yrbs_all_model <- yrbs_all %>%
  
  # Include years of interest
  filter(year %in% c(2001:2021)) %>%
  
  # Exclude children missing age or at least one outcome
  # Note: All other covariates have no missingness due to "Not provided" options
  filter_at(vars(effective_min_wage_l, included_in_study), all_vars(!is.na(.)))

# Rescale weight so mean is 1
yrbs_all_model$weight_2 <- yrbs_all_model$weight/mean(yrbs_all_model$weight, na.rm=T)

# Adolescent characteristics: unweighted
summary(tableby(~ age_2 + sex_2 + race_eth_2 + grade_2,
                yrbs_all_model, digits.pct=1), text=T)

# Adolescent characteristics: weighted
summary(tableby(~ age_2 + sex_2 + race_eth_2 + grade_2,
                yrbs_all_model, digits.pct=1, weights=weight_2), text=T)

# Mental health outcomes: weighted
summary(tableby(~ factor(sad_hopeless), yrbs_all_model, digits.pct=1,
                weights=weight_2, na.action=na.tableby(TRUE)), text=T)
summary(tableby(~ factor(cons_suicide), yrbs_all_model, digits.pct=1,
                weights=weight_2, na.action=na.tableby(TRUE)), text=T)
summary(tableby(~ factor(suicide_att), yrbs_all_model, digits.pct=1,
                weights=weight_2, na.action=na.tableby(TRUE)), text=T)
summary(tableby(~ factor(alcohol), yrbs_all_model, digits.pct=1,
                weights=weight_2, na.action=na.tableby(TRUE)), text=T)
summary(tableby(~ factor(marijuana), yrbs_all_model, digits.pct=1,
                weights=weight_2, na.action=na.tableby(TRUE)), text=T)
summary(tableby(~ factor(fight), yrbs_all_model, digits.pct=1,
                weights=weight_2, na.action=na.tableby(TRUE)), text=T)

# Get missingness by outcome
summary(tableby(~ is.na(sad_hopeless) + is.na(cons_suicide) + is.na(suicide_att) +
                  is.na(alcohol) + is.na(marijuana) + is.na(fight),
                yrbs_all_model, digits.pct=1), text=T)

# State participation by year
write.csv(table(yrbs_all_model$st_abbr.x, yrbs_all_model$year),
          "Exhibits/YRBS state participation.csv")

# Get range of observed minimum wage changes
wage_descriptives <- min_wage_long %>%
  filter(year %in% c(2001:2021), fipsst %in% unique(yrbs_all_model$fipsst)) %>%
  group_by(fipsst) %>%
  summarise(min_wage = min(effective_min_wage_l, na.rm=T),
            max_wage = max(effective_min_wage_l, na.rm=T),
            change   = max_wage - min_wage)

# Range of minimum wages
min(wage_descriptives$min_wage); max(wage_descriptives$max_wage)

# Number of states that never changed
table(wage_descriptives$change > 0)

# Range of changes among states that changed
summary(subset(wage_descriptives, change > 0)$change)

##############################################################################
# Graph: Correlations between outcomes
##############################################################################

# Select desired variables
yrbs_corr <- yrbs_all_model %>% select(sad_hopeless, cons_suicide, suicide_att,
                                         alcohol, marijuana, fight)

# Rename columns
yrbs_corr <- yrbs_corr %>% rename(
  `Sad or hopeless`    = sad_hopeless,
  `Considered suicide` = cons_suicide,
  `Attempted suicide`  = suicide_att,
  `Alcohol`            = alcohol,
  `Marijuana`          = marijuana,
  `Physical fight`     = fight
)

# Run correlation matrix
library(Hmisc); library(corrplot)
corr_matrix <- rcorr(as.matrix(yrbs_corr), type="pearson")

# Create matrices for coefficients, p-values
corr_coef <- corr_matrix$r
corr_pval <- corr_matrix$P

# Define colors for correlation plot
col <- colorRampPalette(c("#4477AA", "#77AADD", "#FFFFFF", "#EE9988", "#BB4444"))

# Plot correlation matrix
corr_plot <- corrplot(corr_coef, method = "color", col = col(200),
                      cl.pos = "n", type = "upper", 
                      
                      # Reorder variables
                      order = "original",
                      
                      # Add coefficient of correlation
                      addCoef.col = "black", number.digits = 2,
                      
                      # Labels color, size, and rotation
                      tl.col = "black", tl.cex = 0.7, number.cex = 0.7, tl.srt = 45,
                      
                      # Hide correlation coefficient on principal diagonal
                      diag = FALSE)
print(corr_plot)

##############################################################################
# Graph: Map of minimum wages
##############################################################################

# Rename states with abbreviations
min_wage_long <- min_wage_long %>% mutate(
  state_map_order = case_when(
    State.or.other.jurisdiction == "Alabama"        ~ "Ala.",
    State.or.other.jurisdiction == "Alaska"         ~ "Alaska",
    State.or.other.jurisdiction == "Arizona"        ~ "Ariz.",
    State.or.other.jurisdiction == "Arkansas"       ~ "Ark.",
    State.or.other.jurisdiction == "California"     ~ "Calif.",
    State.or.other.jurisdiction == "Colorado"       ~ "Colo.",
    State.or.other.jurisdiction == "Connecticut"    ~ "Conn.",
    State.or.other.jurisdiction == "Delaware"       ~ "Del.",
    State.or.other.jurisdiction == "District of Columbia" ~ "D.C.",
    State.or.other.jurisdiction == "Florida"        ~ "Fla.",
    State.or.other.jurisdiction == "Georgia"        ~ "Ga.",
    State.or.other.jurisdiction == "Hawaii"         ~ "Hawaii",
    State.or.other.jurisdiction == "Idaho"          ~ "Idaho",
    State.or.other.jurisdiction == "Illinois"       ~ "Ill.",
    State.or.other.jurisdiction == "Indiana"        ~ "Ind.",
    State.or.other.jurisdiction == "Iowa"           ~ "Iowa",
    State.or.other.jurisdiction == "Kansas"         ~ "Kan.",
    State.or.other.jurisdiction == "Kentucky"       ~ "Ky.",
    State.or.other.jurisdiction == "Louisiana"      ~ "La.",
    State.or.other.jurisdiction == "Maine"          ~ "Maine",
    State.or.other.jurisdiction == "Maryland"       ~ "Md.",
    State.or.other.jurisdiction == "Massachusetts"  ~ "Mass.",
    State.or.other.jurisdiction == "Michigan"       ~ "Mich.",
    State.or.other.jurisdiction == "Minnesota"      ~ "Minn.",
    State.or.other.jurisdiction == "Mississippi"    ~ "Miss.",
    State.or.other.jurisdiction == "Missouri"       ~ "Mo.",
    State.or.other.jurisdiction == "Montana"        ~ "Mont.",
    State.or.other.jurisdiction == "Nebraska"       ~ "Neb.",
    State.or.other.jurisdiction == "Nevada"         ~ "Nev.",
    State.or.other.jurisdiction == "New Hampshire"  ~ "N.H.",
    State.or.other.jurisdiction == "New Jersey"     ~ "N.J.",
    State.or.other.jurisdiction == "New Mexico"     ~ "N.M.",
    State.or.other.jurisdiction == "New York"       ~ "N.Y.",
    State.or.other.jurisdiction == "North Carolina" ~ "N.C.",
    State.or.other.jurisdiction == "North Dakota"   ~ "N.D.",
    State.or.other.jurisdiction == "Ohio"           ~ "Ohio",
    State.or.other.jurisdiction == "Oklahoma"       ~ "Okla.",
    State.or.other.jurisdiction == "Oregon"         ~ "Ore.",
    State.or.other.jurisdiction == "Pennsylvania"   ~ "Pa.",
    State.or.other.jurisdiction == "Rhode Island"   ~ "R.I.",
    State.or.other.jurisdiction == "South Carolina" ~ "S.C.",
    State.or.other.jurisdiction == "South Dakota"   ~ "S.D.",
    State.or.other.jurisdiction == "Tennessee"      ~ "Tenn.",
    State.or.other.jurisdiction == "Texas"          ~ "Texas",
    State.or.other.jurisdiction == "Utah"           ~ "Utah",
    State.or.other.jurisdiction == "Vermont"        ~ "Vt.",
    State.or.other.jurisdiction == "Virginia"       ~ "Va.",
    State.or.other.jurisdiction == "Washington"     ~ "Wash.",
    State.or.other.jurisdiction == "West Virginia"  ~ "W.Va.",
    State.or.other.jurisdiction == "Wisconsin"      ~ "Wis.",
    State.or.other.jurisdiction == "Wyoming"        ~ "Wyo."
  ))

# Generate blanks for empty facets
bl = sapply(1:37, function(n) paste(rep(" ", n), collapse=""))

# Set facet order for states
min_wage_long$state_map_order <-
  factor(min_wage_long$state_map_order,
         levels = c("Alaska", bl[2:10], "Maine", bl[11:19], "Vt.", "N.H.",
                    "Wash.", "Idaho", "Mont.", "N.D.", "Minn.", "Ill.", "Wis.", "Mich.", "N.Y.", "R.I.", "Mass.",
                    "Ore.", "Nev.", "Wyo.", "S.D.", "Iowa", "Ind.", "Ohio", "Pa.", "N.J.", "Conn.", bl[20],
                    "Calif.", "Utah", "Colo.", "Neb.", "Mo.", "Ky.", "W.Va.", "Va.", "Md.", "Del.", bl[21],
                    bl[22], "Ariz.", "N.M.", "Kan.", "Ark.", "Tenn.", "N.C.", "S.C.", "D.C.", bl[23:24],
                    bl[25:27], "Okla.", "La.", "Miss.", "Ala.", "Ga.", bl[28:30],
                    "Hawaii", bl[31:32], "Texas", bl[33:36], "Fla."))

# Subset to study period and states
min_wage_long_sub <- subset(min_wage_long, year %in% c(2001:2022) &
                              fipsst %in% c(1:56))

# Plot minimum wages in shape of U.S.
wage_map_plot <- ggplot(min_wage_long_sub, aes(x=year, y=effective_min_wage_l)) +
  geom_line() +
  facet_wrap(~state_map_order, ncol = 11, drop=F, strip.position="bottom") +
  ylab("Effective minimum wage") +
  xlab("Year (2001–2022)") +
  theme_classic() +
  theme(text = element_text(size = 10, face = "bold"),
        axis.text.x = element_blank(),
        strip.background = element_blank(),
        axis.line  = element_blank(),
        axis.ticks = element_blank()) +
  scale_y_continuous(breaks=c(5.15, 16.10),
                     labels=c("$5.15", "$16.10"))

# Export figure
ggsave(plot=wage_map_plot, file="Exhibits/Minimum wage map plot.pdf",
       width=7, height=5, units='in', dpi=600)

##############################################################################
# Functions for TWFE models
##############################################################################

# Function to extract values from TWFE models.
# Requires: Df for saving coefficients, "lfe" model, title of model.
# Returns: Df of values, ready to pass to "clean_coef_df".
make_coef_df <- function(coef_df, model, TITLE) {
  
  # Get outcome and category labels
  if (model$lhs == "sad_hopeless") {
    outcome  <- "Sad or hopeless in past year"
    category <- "Symptoms"}
  
  if (model$lhs == "cons_suicide") {
    outcome  <- "Considered suicide in past year"
    category <- "Symptoms"}
  
  if (model$lhs == "suicide_att") {
    outcome  <- "Attempted suicide in past year"
    category <- "Symptoms"}
  
  if (model$lhs == "alcohol") {
    outcome  <- "Alcohol use in past month"
    category <- "Substances"}
  
  if (model$lhs == "marijuana") {
    outcome  <- "Marijuana use in past month"
    category <- "Substances"}
  
  if (model$lhs == "fight") {
    outcome  <- "Physical fight in past year"
    category <- "Violence"}
  
  # Add row to coefficient df
  coef_df <- rbind(coef_df, cbind(
    outcome, category, TITLE, model$coef[1], model$cse[1], length(model$residuals)))
  
  # Return coefficient df
  return(coef_df)
}

# Function to clean dataframe of TWFE coefficients.
# Requires: Df with columns specified by "make_coef_df".
# Returns: Cleaned dataframe, ready to use in "print_coef_plot".
clean_coef_df <- function(coef_df) {
  
  # Treat as dataframe
  coef_df <- as.data.frame(coef_df)
  
  # Name columns
  colnames(coef_df) <- c("Outcome", "Category", "Sample", "effect", "se", "n")
  
  # Treat columns as numeric
  coef_df$effect <- as.numeric(coef_df$effect)
  coef_df$se     <- as.numeric(coef_df$se)
  coef_df$n      <- as.numeric(coef_df$n)
  
  # Reorder outcomes
  coef_df$Outcome <- factor(
    coef_df$Outcome, levels=rev(c("Sad or hopeless in past year", "Considered suicide in past year", "Attempted suicide in past year", "Alcohol use in past month", "Marijuana use in past month", "Physical fight in past year")))
  
  # Reorder categories
  coef_df$Category <- factor(
    coef_df$Category, levels=c("Symptoms", "Substances", "Violence"))
  
  # Reorder samples
  coef_df$Sample <- factor(
    coef_df$Sample, levels=rev(c(
      
      # Standard models
      "Adolescents (fully adjusted)",
      "Adolescents (FE only)",
      
      # Sub-population models
      "Black or Hispanic/Latino",
      
      # Lifetime wage models
      "Adolescents (fully adj.), nominal wages",
      "Adolescents (FE only), nominal wages",
      "Adolescents (fully adj.), real wages",
      "Adolescents (FE only), real wages",
      
      # Robustness checks
      "Adolescents (2020 dollars)",
      "Adolescents (lagged wage)",
      
      # Alternate cluster models
      "Adolescents (fully adj., state clusters)",
      "Adolescents (fully adj., nested clusters)",
      "Adolescents (FE only, state clusters)",
      "Adolescents (FE only, nested clusters)"
    )))
  
  # Return df
  return(coef_df)
}

# Function to generate TWFE coefficient plots.
# Requires: Df of coefficients, y title, and y limits.
# Returns: Formatted coefficient plot.
print_coef_plot <- function(coef_df, Y_TITLE, Y_MIN, Y_MAX) {
  
  # Generate coefficient plot
  plot <- ggplot(coef_df, aes(x=Outcome, y=effect, group=Sample, color=Sample)) +
    
    # Null line
    geom_hline(yintercept=0, color="black", linewidth=0.25) +
    
    # Point estimates with distinct shapes
    geom_point(position = position_dodge(width=0.6), size=1, aes(shape=Sample)) +
    scale_shape_manual(values = 1:nlevels(coef_df$Sample)) +
    
    # Confidence intervals
    geom_errorbar(aes(ymin=effect-1.96*se, ymax=effect+1.96*se),
                  position = position_dodge(width=0.6), width=0, linewidth=0.8) +
    geom_errorbar(aes(ymin=effect-2.94*se, ymax=effect+2.94*se),
                  position = position_dodge(width=0.6), width=0.3, alpha=0.5) +
    
    # Titles
    ylab(Y_TITLE) +
    ggtitle("Adolescents (ages 12–18), 2001–2021") +
    
    # Theme modifications
    theme_test() +
    theme(legend.position = "bottom",
          text = element_text(size = 10, face = "bold"),
          plot.title = element_text(hjust=0.5),
          axis.title.y = element_blank(),
          axis.ticks = element_blank(),
          strip.background = element_blank(),
          legend.title = element_blank(),
          panel.grid.major.y = element_blank(),
          panel.grid.major.x = element_line(color="light gray", linewidth=0.5),
          panel.grid.minor.x = element_line(color="light gray", linewidth=0.25)) +
    scale_y_continuous(breaks = seq(-0.1, 0.1, 0.02),
                       minor_breaks = seq(-0.1, 0.1, 0.01),
                       labels = function(x) paste0(x*100," pp")) +
    scale_color_grey(start=0.7, end=0) +
    facet_grid(Category~., scales="free", space="free_y") +
    coord_flip(ylim=c(Y_MIN, Y_MAX)) +
    guides(shape = guide_legend(reverse=T),
           color = guide_legend(reverse=T))
  
  # Return plot
  return(plot)
}

##############################################################################
# Standard TWFE models
##############################################################################

# Sad or hopeless
twfe_min_sad_1 <- felm(sad_hopeless ~ effective_min_wage_l |
                         year + fipsst + birth_year | 0 | fipsst,
                       data    = yrbs_all_model,
                       weights = yrbs_all_model$weight_2)
twfe_min_sad_2 <- felm(sad_hopeless ~ effective_min_wage_l +
                         age_2 + sex_2 + race_eth_2 + grade_2 +
                         elig_1_5 + elig_6_18 + has_eitc + federal_pct + refundable + max_bft_3 |
                         year + fipsst + birth_year | 0 | fipsst,
                       data    = yrbs_all_model,
                       weights = yrbs_all_model$weight_2)

# Considered suicide
twfe_min_con_1 <- felm(cons_suicide ~ effective_min_wage_l |
                         year + fipsst + birth_year | 0 | fipsst,
                       data    = yrbs_all_model,
                       weights = yrbs_all_model$weight_2)
twfe_min_con_2 <- felm(cons_suicide ~ effective_min_wage_l +
                         age_2 + sex_2 + race_eth_2 + grade_2 +
                         elig_1_5 + elig_6_18 + has_eitc + federal_pct + refundable + max_bft_3 |
                         year + fipsst + birth_year | 0 | fipsst,
                       data    = yrbs_all_model,
                       weights = yrbs_all_model$weight_2)

# Attempted suicide
twfe_min_att_1 <- felm(suicide_att ~ effective_min_wage_l |
                         year + fipsst + birth_year | 0 | fipsst,
                       data    = yrbs_all_model,
                       weights = yrbs_all_model$weight_2)
twfe_min_att_2 <- felm(suicide_att ~ effective_min_wage_l +
                         age_2 + sex_2 + race_eth_2 + grade_2 +
                         elig_1_5 + elig_6_18 + has_eitc + federal_pct + refundable + max_bft_3 |
                         year + fipsst + birth_year | 0 | fipsst,
                       data    = yrbs_all_model,
                       weights = yrbs_all_model$weight_2)

# Alcohol use
twfe_min_alc_1 <- felm(alcohol ~ effective_min_wage_l |
                         year + fipsst + birth_year | 0 | fipsst,
                       data    = yrbs_all_model,
                       weights = yrbs_all_model$weight_2)
twfe_min_alc_2 <- felm(alcohol ~ effective_min_wage_l +
                         age_2 + sex_2 + race_eth_2 + grade_2 +
                         elig_1_5 + elig_6_18 + has_eitc + federal_pct + refundable + max_bft_3 |
                         year + fipsst + birth_year | 0 | fipsst,
                       data    = yrbs_all_model,
                       weights = yrbs_all_model$weight_2)

# Marijuana use
twfe_min_mjn_1 <- felm(marijuana ~ effective_min_wage_l |
                         year + fipsst + birth_year | 0 | fipsst,
                       data    = yrbs_all_model,
                       weights = yrbs_all_model$weight_2)
twfe_min_mjn_2 <- felm(marijuana ~ effective_min_wage_l +
                         age_2 + sex_2 + race_eth_2 + grade_2 +
                         elig_1_5 + elig_6_18 + has_eitc + federal_pct + refundable + max_bft_3 |
                         year + fipsst + birth_year | 0 | fipsst,
                       data    = yrbs_all_model,
                       weights = yrbs_all_model$weight_2)

# Physical fight
twfe_min_fgh_1 <- felm(fight ~ effective_min_wage_l |
                         year + fipsst + birth_year | 0 | fipsst,
                       data    = yrbs_all_model,
                       weights = yrbs_all_model$weight_2)
twfe_min_fgh_2 <- felm(fight ~ effective_min_wage_l +
                         age_2 + sex_2 + race_eth_2 + grade_2 +
                         elig_1_5 + elig_6_18 + has_eitc + federal_pct + refundable + max_bft_3 |
                         year + fipsst + birth_year | 0 | fipsst,
                       data    = yrbs_all_model,
                       weights = yrbs_all_model$weight_2)

# Get values from models
twfe_df <- NULL

# Adolescents (FE only)
twfe_df <- make_coef_df(twfe_df, twfe_min_sad_1, "Adolescents (FE only)")
twfe_df <- make_coef_df(twfe_df, twfe_min_con_1, "Adolescents (FE only)")
twfe_df <- make_coef_df(twfe_df, twfe_min_att_1, "Adolescents (FE only)")
twfe_df <- make_coef_df(twfe_df, twfe_min_alc_1, "Adolescents (FE only)")
twfe_df <- make_coef_df(twfe_df, twfe_min_mjn_1, "Adolescents (FE only)")
twfe_df <- make_coef_df(twfe_df, twfe_min_fgh_1, "Adolescents (FE only)")

# Adolescents (fully adjusted)
twfe_df <- make_coef_df(twfe_df, twfe_min_sad_2, "Adolescents (fully adjusted)")
twfe_df <- make_coef_df(twfe_df, twfe_min_con_2, "Adolescents (fully adjusted)")
twfe_df <- make_coef_df(twfe_df, twfe_min_att_2, "Adolescents (fully adjusted)")
twfe_df <- make_coef_df(twfe_df, twfe_min_alc_2, "Adolescents (fully adjusted)")
twfe_df <- make_coef_df(twfe_df, twfe_min_mjn_2, "Adolescents (fully adjusted)")
twfe_df <- make_coef_df(twfe_df, twfe_min_fgh_2, "Adolescents (fully adjusted)")

# Clean dataframe of coefficients
twfe_df <- clean_coef_df(twfe_df)

# Get min. and max. N
min(twfe_df$n); max(twfe_df$n)

# Generate coefficient plot: Standard TWFE (for main text)
plot_twfe_1 <- ggplot(twfe_df, aes(x=Outcome, y=effect, group=Sample, color=Sample)) +
  geom_hline(yintercept=0, color="black", linewidth=0.25) +
  geom_point(position = position_dodge(width=0.6), size=1, aes(shape=Sample)) +
  scale_shape_manual(values = 1:nlevels(twfe_df$Sample)) +
  geom_errorbar(aes(ymin=effect-1.96*se, ymax=effect+1.96*se),
                position = position_dodge(width=0.6), width=0.2) +
  ylab("\nAssociation of $1 increase in minimum wage\nwith rates of mental health outcomes") +
  ggtitle("Adolescents (ages 12–18), 2001–2021") +
  theme_test() +
  theme(legend.position = "bottom",
        text = element_text(size = 10, face = "bold"),
        plot.title = element_text(hjust=0.5),
        axis.title.y = element_blank(),
        axis.ticks = element_blank(),
        strip.background = element_blank(),
        legend.title = element_blank(),
        panel.grid.major.y = element_blank(),
        panel.grid.major.x = element_line(color="light gray", linewidth=0.5),
        panel.grid.minor.x = element_line(color="light gray", linewidth=0.25)) +
  scale_y_continuous(breaks = seq(-0.1, 0.1, 0.02),
                     minor_breaks = seq(-0.1, 0.1, 0.01),
                     labels = function(x) paste0(x*100," pp")) +
  scale_color_grey(start=0.7, end=0) +
  facet_grid(Category~., scales="free", space="free_y") +
  coord_flip(ylim=c(-0.045, 0.045)) +
  guides(shape = guide_legend(reverse=T),
         color = guide_legend(reverse=T))

# Export figure
ggsave(plot=plot_twfe_1, file="Exhibits/YRBS coefficient plot, TWFE 1.pdf",
       width=5.75, height=4.5, units='in', dpi=600)

# Generate coefficient plot: Standard TWFE (for appendix)
plot_twfe_2 <- print_coef_plot(
  twfe_df,
  Y_TITLE    = "\nAssociation of $1 increase in minimum wage\nwith rates of mental health outcomes",
  Y_MIN      = -0.045,
  Y_MAX      =  0.045
  )

# Export figure
ggsave(plot=plot_twfe_2, file="Exhibits/YRBS coefficient plot, TWFE 2.pdf",
       width=5.75, height=4.5, units='in', dpi=600)

# Variable names for table
var_names <- 
  c("effective_min_wage_l" = "$1 increase in minimum wage",
    "age_213 years old" = "13 years old",
    "age_214 years old" = "14 years old",
    "age_215 years old" = "15 years old",
    "age_216 years old" = "16 years old",
    "age_217 years old" = "17 years old",
    "age_218 years old or older" = "18 years old or older",
    "sex_2Female" = "Female",
    "sex_2Not provided" = "Sex not provided",
    "race_eth_2American Indian or Alaska Native" =
      "American Indian or Alaska Native",
    "race_eth_2Asian, Native Hawaiian, or Pacific Islander" =
      "Asian, Native Hawaiian, or Pacific Islander",
    "race_eth_2Black or African American" = "Black or African American",
    "race_eth_2Hispanic/Latino" = "Hispanic/Latino",
    "race_eth_2White" = "White",
    "race_eth_2Multiple races, non-Hispanic" = "Multiple races, non-Hispanic",
    "race_eth_2Not provided" = "Race/ethnicity not provided",
    "grade_210th grade" = "10th grade",
    "grade_211th grade" = "11th grade",
    "grade_212th grade" = "12th grade",
    "grade_2Not provided" = "Grade not provided",
    "elig_1_5"    = "Medicare income eligibility limit, age 1-5",
    "elig_6_18"   = "Medicare income eligibility limit, age 6-18",
    "has_eitc"    = "State has EITC",
    "federal_pct" = "State EITC as percent of federal",
    "refundable"  = "State EITC is refundable",
    "max_bft_3"   = "Maximum TANF benefit for family of 3")

# Compile results into tables
modelsummary(list("FE only"    = twfe_min_sad_1,
                  "Fully adj." = twfe_min_sad_2,
                  "FE only"    = twfe_min_con_1,
                  "Fully adj." = twfe_min_con_2,
                  "FE only"    = twfe_min_att_1,
                  "Fully adj." = twfe_min_att_2,
                  "FE only"    = twfe_min_alc_1,
                  "Fully adj." = twfe_min_alc_2,
                  "FE only"    = twfe_min_mjn_1,
                  "Fully adj." = twfe_min_mjn_2,
                  "FE only"    = twfe_min_fgh_1,
                  "Fully adj." = twfe_min_fgh_2),
             gof_omit    = "Log*|AIC|BIC|F|RMSE|Std",
             
             # Options for viewing abbreviated table in R
             coef_omit   = "^(?!effective_min_wage_l)",
             statistic   = c("[{conf.low} to {conf.high}]", "P={p.value}"),
             fmt         = fmt_statistic("estimate"=4, "conf.low"=4, "conf.high"=4, "p.value"=2),
             
             # Options for exporting full table as .csv
             # coef_rename = var_names,
             # statistic   = c("(({std.error}))", "P={p.value}"),
             # fmt         = fmt_statistic("estimate"=4, "std.error"=4, "p.value"=3),
             # output      = "YRBS table, standard TWFE.csv"
) %>%
  add_header_above(c(" "=1, "Sad or hopeless"=2, "Considered suicide"=2,
                     "Attempted suicide"=2, "Alcohol"=2, "Marijuana"=2, "Physical fights"=2))

##############################################################################
# TWFE robustness check: Subgroup
##############################################################################

# Generate subgroup of interest
yrbs_all_model_r <- subset(yrbs_all_model, race7 %in% c(3:4)) # Black or Hispanic/Latino

# Sad or hopeless
twfe_min_sad_r <- felm(sad_hopeless ~ effective_min_wage_l +
                         age_2 + sex_2 + race_eth_2 + grade_2 +
                         elig_1_5 + elig_6_18 + has_eitc + federal_pct + refundable + max_bft_3 |
                         year + fipsst + birth_year | 0 | fipsst,
                       data    = yrbs_all_model_r,
                       weights = yrbs_all_model_r$weight_2)

# Considered suicide
twfe_min_con_r <- felm(cons_suicide ~ effective_min_wage_l +
                         age_2 + sex_2 + race_eth_2 + grade_2 +
                         elig_1_5 + elig_6_18 + has_eitc + federal_pct + refundable + max_bft_3 |
                         year + fipsst + birth_year | 0 | fipsst,
                       data    = yrbs_all_model_r,
                       weights = yrbs_all_model_r$weight_2)

# Attempted suicide
twfe_min_att_r <- felm(suicide_att ~ effective_min_wage_l +
                         age_2 + sex_2 + race_eth_2 + grade_2 +
                         elig_1_5 + elig_6_18 + has_eitc + federal_pct + refundable + max_bft_3 |
                         year + fipsst + birth_year | 0 | fipsst,
                       data    = yrbs_all_model_r,
                       weights = yrbs_all_model_r$weight_2)

# Alcohol use
twfe_min_alc_r <- felm(alcohol ~ effective_min_wage_l +
                         age_2 + sex_2 + race_eth_2 + grade_2 +
                         elig_1_5 + elig_6_18 + has_eitc + federal_pct + refundable + max_bft_3 |
                         year + fipsst + birth_year | 0 | fipsst,
                       data    = yrbs_all_model_r,
                       weights = yrbs_all_model_r$weight_2)

# Marijuana use
twfe_min_mjn_r <- felm(marijuana ~ effective_min_wage_l +
                         age_2 + sex_2 + race_eth_2 + grade_2 +
                         elig_1_5 + elig_6_18 + has_eitc + federal_pct + refundable + max_bft_3 |
                         year + fipsst + birth_year | 0 | fipsst,
                       data    = yrbs_all_model_r,
                       weights = yrbs_all_model_r$weight_2)

# Physical fight
twfe_min_fgh_r <- felm(fight ~ effective_min_wage_l +
                         age_2 + sex_2 + race_eth_2 + grade_2 +
                         elig_1_5 + elig_6_18 + has_eitc + federal_pct + refundable + max_bft_3 |
                         year + fipsst + birth_year | 0 | fipsst,
                       data    = yrbs_all_model_r,
                       weights = yrbs_all_model_r$weight_2)

# Get values from models
sub_df <- NULL

# Black or Hispanic/Latino
sub_df <- make_coef_df(sub_df, twfe_min_sad_r, "Black or Hispanic/Latino")
sub_df <- make_coef_df(sub_df, twfe_min_con_r, "Black or Hispanic/Latino")
sub_df <- make_coef_df(sub_df, twfe_min_att_r, "Black or Hispanic/Latino")
sub_df <- make_coef_df(sub_df, twfe_min_alc_r, "Black or Hispanic/Latino")
sub_df <- make_coef_df(sub_df, twfe_min_mjn_r, "Black or Hispanic/Latino")
sub_df <- make_coef_df(sub_df, twfe_min_fgh_r, "Black or Hispanic/Latino")

# Clean dataframe of coefficients
sub_df <- clean_coef_df(sub_df)

# Get min. and max. N
# Before adding standard TWFE models
min(sub_df$n); max(sub_df$n)

# Add standard TWFE models for comparison
sub_df <- rbind(sub_df, twfe_df %>% filter(Sample == "Adolescents (fully adjusted)"))

# Generate coefficient plot: Subgroup
plot_sub <- print_coef_plot(
  sub_df,
  Y_TITLE    = "\nAssociation of $1 increase in minimum wage\nwith rates of mental health outcomes",
  Y_MIN      = -0.045,
  Y_MAX      =  0.045
  )

# Adjust legend
plot_sub <- plot_sub +
  guides(shape = guide_legend(reverse=T, nrow=2),
         color = guide_legend(reverse=T, nrow=2))

# Export figure
ggsave(plot=plot_sub, file="Exhibits/YRBS coefficient plot, subgroup.pdf",
       width=5.75, height=4.5, units='in', dpi=600)

##############################################################################
# TWFE robustness check: Alternate specifications
##############################################################################

# Sad or hopeless
twfe_min_sad_x <- felm(sad_hopeless ~ inflation_min_wage_l +
                         age_2 + sex_2 + race_eth_2 + grade_2 +
                         elig_1_5 + elig_6_18 + has_eitc + federal_pct + refundable + max_bft_3 |
                         year + fipsst + birth_year | 0 | fipsst,
                       data    = yrbs_all_model,
                       weights = yrbs_all_model$weight_2)
twfe_min_sad_y <- felm(sad_hopeless ~ lag_by_1 +
                         age_2 + sex_2 + race_eth_2 + grade_2 +
                         elig_1_5 + elig_6_18 + has_eitc + federal_pct + refundable + max_bft_3 |
                         year + fipsst + birth_year | 0 | fipsst,
                       data    = yrbs_all_model,
                       weights = yrbs_all_model$weight_2)

# Considered suicide
twfe_min_con_x <- felm(cons_suicide ~ inflation_min_wage_l +
                         age_2 + sex_2 + race_eth_2 + grade_2 +
                         elig_1_5 + elig_6_18 + has_eitc + federal_pct + refundable + max_bft_3 |
                         year + fipsst + birth_year | 0 | fipsst,
                       data    = yrbs_all_model,
                       weights = yrbs_all_model$weight_2)
twfe_min_con_y <- felm(cons_suicide ~ lag_by_1 +
                         age_2 + sex_2 + race_eth_2 + grade_2 +
                         elig_1_5 + elig_6_18 + has_eitc + federal_pct + refundable + max_bft_3 |
                         year + fipsst + birth_year | 0 | fipsst,
                       data    = yrbs_all_model,
                       weights = yrbs_all_model$weight_2)

# Attempted suicide
twfe_min_att_x <- felm(suicide_att ~ inflation_min_wage_l +
                         age_2 + sex_2 + race_eth_2 + grade_2 +
                         elig_1_5 + elig_6_18 + has_eitc + federal_pct + refundable + max_bft_3 |
                         year + fipsst + birth_year | 0 | fipsst,
                       data    = yrbs_all_model,
                       weights = yrbs_all_model$weight_2)
twfe_min_att_y <- felm(suicide_att ~ lag_by_1 +
                         age_2 + sex_2 + race_eth_2 + grade_2 +
                         elig_1_5 + elig_6_18 + has_eitc + federal_pct + refundable + max_bft_3 |
                         year + fipsst + birth_year | 0 | fipsst,
                       data    = yrbs_all_model,
                       weights = yrbs_all_model$weight_2)

# Alcohol use
twfe_min_alc_x <- felm(alcohol ~ inflation_min_wage_l +
                         age_2 + sex_2 + race_eth_2 + grade_2 +
                         elig_1_5 + elig_6_18 + has_eitc + federal_pct + refundable + max_bft_3 |
                         year + fipsst + birth_year | 0 | fipsst,
                       data    = yrbs_all_model,
                       weights = yrbs_all_model$weight_2)
twfe_min_alc_y <- felm(alcohol ~ lag_by_1 +
                         age_2 + sex_2 + race_eth_2 + grade_2 +
                         elig_1_5 + elig_6_18 + has_eitc + federal_pct + refundable + max_bft_3 |
                         year + fipsst + birth_year | 0 | fipsst,
                       data    = yrbs_all_model,
                       weights = yrbs_all_model$weight_2)

# Marijuana use
twfe_min_mjn_x <- felm(marijuana ~ inflation_min_wage_l +
                         age_2 + sex_2 + race_eth_2 + grade_2 +
                         elig_1_5 + elig_6_18 + has_eitc + federal_pct + refundable + max_bft_3 |
                         year + fipsst + birth_year | 0 | fipsst,
                       data    = yrbs_all_model,
                       weights = yrbs_all_model$weight_2)
twfe_min_mjn_y <- felm(marijuana ~ lag_by_1 +
                         age_2 + sex_2 + race_eth_2 + grade_2 +
                         elig_1_5 + elig_6_18 + has_eitc + federal_pct + refundable + max_bft_3 |
                         year + fipsst + birth_year | 0 | fipsst,
                       data    = yrbs_all_model,
                       weights = yrbs_all_model$weight_2)

# Physical fight
twfe_min_fgh_x <- felm(fight ~ inflation_min_wage_l +
                         age_2 + sex_2 + race_eth_2 + grade_2 +
                         elig_1_5 + elig_6_18 + has_eitc + federal_pct + refundable + max_bft_3 |
                         year + fipsst + birth_year | 0 | fipsst,
                       data    = yrbs_all_model,
                       weights = yrbs_all_model$weight_2)
twfe_min_fgh_y <- felm(fight ~ lag_by_1 +
                         age_2 + sex_2 + race_eth_2 + grade_2 +
                         elig_1_5 + elig_6_18 + has_eitc + federal_pct + refundable + max_bft_3 |
                         year + fipsst + birth_year | 0 | fipsst,
                       data    = yrbs_all_model,
                       weights = yrbs_all_model$weight_2)

# Get values from models
rob_df <- NULL

# Adolescents (2020 dollars)
rob_df <- make_coef_df(rob_df, twfe_min_sad_x, "Adolescents (2020 dollars)")
rob_df <- make_coef_df(rob_df, twfe_min_con_x, "Adolescents (2020 dollars)")
rob_df <- make_coef_df(rob_df, twfe_min_att_x, "Adolescents (2020 dollars)")
rob_df <- make_coef_df(rob_df, twfe_min_alc_x, "Adolescents (2020 dollars)")
rob_df <- make_coef_df(rob_df, twfe_min_mjn_x, "Adolescents (2020 dollars)")
rob_df <- make_coef_df(rob_df, twfe_min_fgh_x, "Adolescents (2020 dollars)")

# Adolescents (lagged wage)
rob_df <- make_coef_df(rob_df, twfe_min_sad_y, "Adolescents (lagged wage)")
rob_df <- make_coef_df(rob_df, twfe_min_con_y, "Adolescents (lagged wage)")
rob_df <- make_coef_df(rob_df, twfe_min_att_y, "Adolescents (lagged wage)")
rob_df <- make_coef_df(rob_df, twfe_min_alc_y, "Adolescents (lagged wage)")
rob_df <- make_coef_df(rob_df, twfe_min_mjn_y, "Adolescents (lagged wage)")
rob_df <- make_coef_df(rob_df, twfe_min_fgh_y, "Adolescents (lagged wage)")

# Clean dataframe of coefficients
rob_df <- clean_coef_df(rob_df)

# Add standard TWFE models for comparison
rob_df <- rbind(rob_df, twfe_df %>% filter(Sample == "Adolescents (fully adjusted)"))

# Get min. and max. N
min(rob_df$n); max(rob_df$n)

# Generate coefficient plot: Alternate specifications
plot_rob <- print_coef_plot(
  rob_df,
  Y_TITLE    = "\nAssociation of $1 increase in minimum wage\nwith rates of mental health outcomes",
  Y_MIN      = -0.045,
  Y_MAX      =  0.045
  )

# Adjust legend
plot_rob <- plot_rob +
  guides(shape = guide_legend(reverse=T, nrow=3),
         color = guide_legend(reverse=T, nrow=3))

# Export figure
ggsave(plot=plot_rob, file="Exhibits/YRBS coefficient plot, robustness.pdf",
       width=5.75, height=5, units='in', dpi=600)

##############################################################################
# TWFE robustness check: Logistic regression
##############################################################################

# Treat fixed effects as factors
yrbs_all_model$year       <- as.factor(yrbs_all_model$year)
yrbs_all_model$fipsst     <- as.factor(yrbs_all_model$fipsst)
yrbs_all_model$birth_year <- as.factor(yrbs_all_model$birth_year)

# Define complex sampling design
# Use state clusters and apply survey weights
design_yrbs <- svydesign(id=~fipsst, weights=~weight_2, data=yrbs_all_model)

# Sad or hopeless
log_min_sad_1 <- svyglm(sad_hopeless ~ effective_min_wage_l +
                          year + fipsst + birth_year,
                        design = design_yrbs, family="quasibinomial")
log_min_sad_2 <- svyglm(sad_hopeless ~ effective_min_wage_l +
                          age_2 + sex_2 + race_eth_2 + grade_2 +
                          elig_1_5 + elig_6_18 + has_eitc + federal_pct + refundable + max_bft_3 +
                          year + fipsst + birth_year,
                        design = design_yrbs, family="quasibinomial")

# Considered suicide
log_min_con_1 <- svyglm(cons_suicide ~ effective_min_wage_l +
                          year + fipsst + birth_year,
                        design = design_yrbs, family="quasibinomial")
log_min_con_2 <- svyglm(cons_suicide ~ effective_min_wage_l +
                          age_2 + sex_2 + race_eth_2 + grade_2 +
                          elig_1_5 + elig_6_18 + has_eitc + federal_pct + refundable + max_bft_3 +
                          year + fipsst + birth_year,
                        design = design_yrbs, family="quasibinomial")

# Attempted suicide
log_min_att_1 <- svyglm(suicide_att ~ effective_min_wage_l +
                          year + fipsst + birth_year,
                        design = design_yrbs, family="quasibinomial")
log_min_att_2 <- svyglm(suicide_att ~ effective_min_wage_l +
                          age_2 + sex_2 + race_eth_2 + grade_2 +
                          elig_1_5 + elig_6_18 + has_eitc + federal_pct + refundable + max_bft_3 +
                          year + fipsst + birth_year,
                        design = design_yrbs, family="quasibinomial")

# Alcohol use
log_min_alc_1 <- svyglm(alcohol ~ effective_min_wage_l +
                          year + fipsst + birth_year,
                        design = design_yrbs, family="quasibinomial")
log_min_alc_2 <- svyglm(alcohol ~ effective_min_wage_l +
                          age_2 + sex_2 + race_eth_2 + grade_2 +
                          elig_1_5 + elig_6_18 + has_eitc + federal_pct + refundable + max_bft_3 +
                          year + fipsst + birth_year,
                        design = design_yrbs, family="quasibinomial")

# Marijuana use
log_min_mjn_1 <- svyglm(marijuana ~ effective_min_wage_l +
                          year + fipsst + birth_year,
                        design = design_yrbs, family="quasibinomial")
log_min_mjn_2 <- svyglm(marijuana ~ effective_min_wage_l +
                          age_2 + sex_2 + race_eth_2 + grade_2 +
                          elig_1_5 + elig_6_18 + has_eitc + federal_pct + refundable + max_bft_3 +
                          year + fipsst + birth_year,
                        design = design_yrbs, family="quasibinomial")

# Physical fight
log_min_fgh_1 <- svyglm(fight ~ effective_min_wage_l +
                          year + fipsst + birth_year,
                        design = design_yrbs, family="quasibinomial")
log_min_fgh_2 <- svyglm(fight ~ effective_min_wage_l +
                          age_2 + sex_2 + race_eth_2 + grade_2 +
                          elig_1_5 + elig_6_18 + has_eitc + federal_pct + refundable + max_bft_3 +
                          year + fipsst + birth_year,
                        design = design_yrbs, family="quasibinomial")

# Get values from models
log_df <- as.data.frame(rbind(
  # Sad or hopeless
  cbind("Sad or hopeless in past year", "Symptoms", "Adolescents (FE only)",
        coef(log_min_sad_1)[2],
        SE(log_min_sad_1)[2],
        length(log_min_sad_1$residuals)),
  
  cbind("Sad or hopeless in past year", "Symptoms", "Adolescents (fully adjusted)",
        coef(log_min_sad_2)[2],
        SE(log_min_sad_2)[2],
        length(log_min_sad_2$residuals)),
  
  # Considered suicide
  cbind("Considered suicide in past year", "Symptoms", "Adolescents (FE only)",
        coef(log_min_con_1)[2],
        SE(log_min_con_1)[2],
        length(log_min_con_1$residuals)),
  
  cbind("Considered suicide in past year", "Symptoms", "Adolescents (fully adjusted)",
        coef(log_min_con_2)[2],
        SE(log_min_con_2)[2],
        length(log_min_con_2$residuals)),
  
  # Attempted suicide
  cbind("Attempted suicide in past year", "Symptoms", "Adolescents (FE only)",
        coef(log_min_att_1)[2],
        SE(log_min_att_1)[2],
        length(log_min_att_1$residuals)),
  
  cbind("Attempted suicide in past year", "Symptoms", "Adolescents (fully adjusted)",
        coef(log_min_att_2)[2],
        SE(log_min_att_2)[2],
        length(log_min_att_2$residuals)),
  
  # Alcohol use
  cbind("Alcohol use in past month", "Substances", "Adolescents (FE only)",
        coef(log_min_alc_1)[2],
        SE(log_min_alc_1)[2],
        length(log_min_alc_1$residuals)),
  
  cbind("Alcohol use in past month", "Substances", "Adolescents (fully adjusted)",
        coef(log_min_alc_2)[2],
        SE(log_min_alc_2)[2],
        length(log_min_alc_2$residuals)),
  
  # Marijuana use
  cbind("Marijuana use in past month", "Substances", "Adolescents (FE only)",
        coef(log_min_mjn_1)[2],
        SE(log_min_mjn_1)[2],
        length(log_min_mjn_1$residuals)),
  
  cbind("Marijuana use in past month", "Substances", "Adolescents (fully adjusted)",
        coef(log_min_mjn_2)[2],
        SE(log_min_mjn_2)[2],
        length(log_min_mjn_2$residuals)),
  
  # Physical fight
  cbind("Physical fight in past year", "Violence", "Adolescents (FE only)",
        coef(log_min_fgh_1)[2],
        SE(log_min_fgh_1)[2],
        length(log_min_fgh_1$residuals)),
  
  cbind("Physical fight in past year", "Violence", "Adolescents (fully adjusted)",
        coef(log_min_fgh_2)[2],
        SE(log_min_fgh_2)[2],
        length(log_min_fgh_2$residuals))
))

# Name columns
colnames(log_df) <- c("Outcome", "Category", "Sample", "log_odds", "se", "n")

# Reorder outcomes
log_df$Outcome <- factor(
  log_df$Outcome, levels=c("Sad or hopeless in past year", "Considered suicide in past year", "Attempted suicide in past year", "Alcohol use in past month", "Marijuana use in past month", "Physical fight in past year"))

# Reorder categories
log_df$Category <- factor(
  log_df$Category, levels=c("Symptoms", "Substances", "Violence"))

# Reorder samples
log_df$Sample <- factor(log_df$Sample, levels=rev(c("Adolescents (fully adjusted)",
                                                    "Adolescents (FE only)")))

# Treat columns as numeric
log_df$log_odds <- as.numeric(log_df$log_odds)
log_df$se       <- as.numeric(log_df$se)
log_df$n        <- as.numeric(log_df$n)

# Get min. and max. N
min(log_df$n); max(log_df$n)

# Generate coefficient plot: Logistic
plot_log <- ggplot(log_df, aes(x=Outcome, y=exp(log_odds), group=Sample, color=Sample)) +
  geom_hline(yintercept=0, color="black", linewidth=0.25) +
  geom_point(position = position_dodge(width=0.6), size=1, aes(shape=Sample)) +
  scale_shape_manual(values = 1:nlevels(log_df$Sample)) +
  geom_errorbar(aes(ymin=exp(log_odds - 1.96*se), ymax=exp(log_odds + 1.96*se)),
                position = position_dodge(width=0.6), width=0, linewidth=0.8) +
  geom_errorbar(aes(ymin=exp(log_odds - 2.94*se), ymax=exp(log_odds + 2.94*se)),
                position = position_dodge(width=0.6), width=0.3, alpha=0.5) +
  ylab("\nAssociation of $1 increase in minimum wage\nwith odds of mental health outcomes") +
  ggtitle("Adolescents (ages 12–18), 2001–2021") +
  theme_test() +
  theme(legend.position = "bottom",
        text = element_text(size = 10, face = "bold"),
        plot.title = element_text(hjust=0.5),
        axis.title.y = element_blank(),
        axis.ticks = element_blank(),
        strip.background = element_blank(),
        legend.title = element_blank(),
        panel.grid.major.y = element_blank(),
        panel.grid.major.x = element_line(color="light gray", linewidth=0.5),
        panel.grid.minor.x = element_line(color="light gray", linewidth=0.25)) +
  scale_y_continuous(trans = "log10", limits=c(0.25,4),
                     labels=c("0.25", "0.5", "1\nOdds ratio", "2", "4"),
                     breaks=c(0.25, 0.5, 1, 2, 4),
                     minor_breaks = seq(0.25, 10, 0.25)) +
  scale_color_grey(start=0.7, end=0) +
  facet_grid(Category~., scales="free", space="free_y") +
  coord_flip() +
  guides(shape = guide_legend(reverse=T),
         color = guide_legend(reverse=T))

# Export figure
ggsave(plot=plot_log, file="Exhibits/YRBS coefficient plot, logistic.pdf",
       width=5.75, height=4.5, units='in', dpi=600)

##############################################################################
# TWFE robustness check: Level-log models
##############################################################################

# Sad or hopeless
twfe_min_sad_ln_1 <- felm(sad_hopeless ~ log(effective_min_wage_l) |
                            year + fipsst + birth_year | 0 | fipsst,
                          data    = yrbs_all_model,
                          weights = yrbs_all_model$weight_2)
twfe_min_sad_ln_2 <- felm(sad_hopeless ~ log(effective_min_wage_l) +
                            age_2 + sex_2 + race_eth_2 + grade_2 +
                            elig_1_5 + elig_6_18 + has_eitc + federal_pct + refundable + max_bft_3 |
                            year + fipsst + birth_year | 0 | fipsst,
                          data    = yrbs_all_model,
                          weights = yrbs_all_model$weight_2)

# Considered suicide
twfe_min_con_ln_1 <- felm(cons_suicide ~ log(effective_min_wage_l) |
                            year + fipsst + birth_year | 0 | fipsst,
                          data    = yrbs_all_model,
                          weights = yrbs_all_model$weight_2)
twfe_min_con_ln_2 <- felm(cons_suicide ~ log(effective_min_wage_l) +
                            age_2 + sex_2 + race_eth_2 + grade_2 +
                            elig_1_5 + elig_6_18 + has_eitc + federal_pct + refundable + max_bft_3 |
                            year + fipsst + birth_year | 0 | fipsst,
                          data    = yrbs_all_model,
                          weights = yrbs_all_model$weight_2)

# Attempted suicide
twfe_min_att_ln_1 <- felm(suicide_att ~ log(effective_min_wage_l) |
                            year + fipsst + birth_year | 0 | fipsst,
                          data    = yrbs_all_model,
                          weights = yrbs_all_model$weight_2)
twfe_min_att_ln_2 <- felm(suicide_att ~ log(effective_min_wage_l) +
                            age_2 + sex_2 + race_eth_2 + grade_2 +
                            elig_1_5 + elig_6_18 + has_eitc + federal_pct + refundable + max_bft_3 |
                            year + fipsst + birth_year | 0 | fipsst,
                          data    = yrbs_all_model,
                          weights = yrbs_all_model$weight_2)

# Alcohol use
twfe_min_alc_ln_1 <- felm(alcohol ~ log(effective_min_wage_l) |
                            year + fipsst + birth_year | 0 | fipsst,
                          data    = yrbs_all_model,
                          weights = yrbs_all_model$weight_2)
twfe_min_alc_ln_2 <- felm(alcohol ~ log(effective_min_wage_l) +
                            age_2 + sex_2 + race_eth_2 + grade_2 +
                            elig_1_5 + elig_6_18 + has_eitc + federal_pct + refundable + max_bft_3 |
                            year + fipsst + birth_year | 0 | fipsst,
                          data    = yrbs_all_model,
                          weights = yrbs_all_model$weight_2)

# Marijuana use
twfe_min_mjn_ln_1 <- felm(marijuana ~ log(effective_min_wage_l) |
                            year + fipsst + birth_year | 0 | fipsst,
                          data    = yrbs_all_model,
                          weights = yrbs_all_model$weight_2)
twfe_min_mjn_ln_2 <- felm(marijuana ~ log(effective_min_wage_l) +
                            age_2 + sex_2 + race_eth_2 + grade_2 +
                            elig_1_5 + elig_6_18 + has_eitc + federal_pct + refundable + max_bft_3 |
                            year + fipsst + birth_year | 0 | fipsst,
                          data    = yrbs_all_model,
                          weights = yrbs_all_model$weight_2)

# Physical fight
twfe_min_fgh_ln_1 <- felm(fight ~ log(effective_min_wage_l) |
                            year + fipsst + birth_year | 0 | fipsst,
                          data    = yrbs_all_model,
                          weights = yrbs_all_model$weight_2)
twfe_min_fgh_ln_2 <- felm(fight ~ log(effective_min_wage_l) +
                            age_2 + sex_2 + race_eth_2 + grade_2 +
                            elig_1_5 + elig_6_18 + has_eitc + federal_pct + refundable + max_bft_3 |
                            year + fipsst + birth_year | 0 | fipsst,
                          data    = yrbs_all_model,
                          weights = yrbs_all_model$weight_2)

# Get values from models
level_log_df <- NULL

# Adolescents (FE only)
level_log_df <- make_coef_df(level_log_df, twfe_min_sad_ln_1, "Adolescents (FE only)")
level_log_df <- make_coef_df(level_log_df, twfe_min_con_ln_1, "Adolescents (FE only)")
level_log_df <- make_coef_df(level_log_df, twfe_min_att_ln_1, "Adolescents (FE only)")
level_log_df <- make_coef_df(level_log_df, twfe_min_alc_ln_1, "Adolescents (FE only)")
level_log_df <- make_coef_df(level_log_df, twfe_min_mjn_ln_1, "Adolescents (FE only)")
level_log_df <- make_coef_df(level_log_df, twfe_min_fgh_ln_1, "Adolescents (FE only)")

# Adolescents (fully adjusted)
level_log_df <- make_coef_df(level_log_df, twfe_min_sad_ln_2, "Adolescents (fully adjusted)")
level_log_df <- make_coef_df(level_log_df, twfe_min_con_ln_2, "Adolescents (fully adjusted)")
level_log_df <- make_coef_df(level_log_df, twfe_min_att_ln_2, "Adolescents (fully adjusted)")
level_log_df <- make_coef_df(level_log_df, twfe_min_alc_ln_2, "Adolescents (fully adjusted)")
level_log_df <- make_coef_df(level_log_df, twfe_min_mjn_ln_2, "Adolescents (fully adjusted)")
level_log_df <- make_coef_df(level_log_df, twfe_min_fgh_ln_2, "Adolescents (fully adjusted)")

# Clean dataframe of coefficients
level_log_df <- clean_coef_df(level_log_df)

# Rescale coefficients for 10% increase in minimum wage
level_log_df$effect <- level_log_df$effect*10*0.01
level_log_df$se     <- level_log_df$se*10*0.01

# Get min. and max. N
min(level_log_df$n); max(level_log_df$n)

# Generate coefficient plot: Alternate specifications
plot_level_log <- print_coef_plot(
  level_log_df,
  Y_TITLE    = "\nAssociation of 10% increase in minimum wage\nwith rates of mental health outcomes",
  Y_MIN      = -0.045,
  Y_MAX      =  0.045
)

# Export figure
ggsave(plot=plot_level_log, file="Exhibits/YRBS coefficient plot, level-log.pdf",
       width=5.75, height=4.5, units='in', dpi=600)

##############################################################################
# TWFE robustness check: Nested clusters
##############################################################################

# Sad or hopeless
twfe_min_sad_1c <- felm(sad_hopeless ~ effective_min_wage_l |
                          year + fipsst + birth_year | 0 | cluster,
                        data    = yrbs_all_model,
                        weights = yrbs_all_model$weight_2)
twfe_min_sad_2c <- felm(sad_hopeless ~ effective_min_wage_l +
                          age_2 + sex_2 + race_eth_2 + grade_2 +
                          elig_1_5 + elig_6_18 + has_eitc + federal_pct + refundable + max_bft_3 |
                          year + fipsst + birth_year | 0 | cluster,
                        data    = yrbs_all_model,
                        weights = yrbs_all_model$weight_2)

# Considered suicide
twfe_min_con_1c <- felm(cons_suicide ~ effective_min_wage_l |
                          year + fipsst + birth_year | 0 | cluster,
                        data    = yrbs_all_model,
                        weights = yrbs_all_model$weight_2)
twfe_min_con_2c <- felm(cons_suicide ~ effective_min_wage_l +
                          age_2 + sex_2 + race_eth_2 + grade_2 +
                          elig_1_5 + elig_6_18 + has_eitc + federal_pct + refundable + max_bft_3 |
                          year + fipsst + birth_year | 0 | cluster,
                        data    = yrbs_all_model,
                        weights = yrbs_all_model$weight_2)

# Attempted suicide
twfe_min_att_1c <- felm(suicide_att ~ effective_min_wage_l |
                          year + fipsst + birth_year | 0 | cluster,
                        data    = yrbs_all_model,
                        weights = yrbs_all_model$weight_2)
twfe_min_att_2c <- felm(suicide_att ~ effective_min_wage_l +
                          age_2 + sex_2 + race_eth_2 + grade_2 +
                          elig_1_5 + elig_6_18 + has_eitc + federal_pct + refundable + max_bft_3 |
                          year + fipsst + birth_year | 0 | cluster,
                        data    = yrbs_all_model,
                        weights = yrbs_all_model$weight_2)

# Alcohol use
twfe_min_alc_1c <- felm(alcohol ~ effective_min_wage_l |
                          year + fipsst + birth_year | 0 | cluster,
                        data    = yrbs_all_model,
                        weights = yrbs_all_model$weight_2)
twfe_min_alc_2c <- felm(alcohol ~ effective_min_wage_l +
                          age_2 + sex_2 + race_eth_2 + grade_2 +
                          elig_1_5 + elig_6_18 + has_eitc + federal_pct + refundable + max_bft_3 |
                          year + fipsst + birth_year | 0 | cluster,
                        data    = yrbs_all_model,
                        weights = yrbs_all_model$weight_2)

# Marijuana use
twfe_min_mjn_1c <- felm(marijuana ~ effective_min_wage_l |
                          year + fipsst + birth_year | 0 | cluster,
                        data    = yrbs_all_model,
                        weights = yrbs_all_model$weight_2)
twfe_min_mjn_2c <- felm(marijuana ~ effective_min_wage_l +
                          age_2 + sex_2 + race_eth_2 + grade_2 +
                          elig_1_5 + elig_6_18 + has_eitc + federal_pct + refundable + max_bft_3 |
                          year + fipsst + birth_year | 0 | cluster,
                        data    = yrbs_all_model,
                        weights = yrbs_all_model$weight_2)

# Physical fight
twfe_min_fgh_1c <- felm(fight ~ effective_min_wage_l |
                          year + fipsst + birth_year | 0 | cluster,
                        data    = yrbs_all_model,
                        weights = yrbs_all_model$weight_2)
twfe_min_fgh_2c <- felm(fight ~ effective_min_wage_l +
                          age_2 + sex_2 + race_eth_2 + grade_2 +
                          elig_1_5 + elig_6_18 + has_eitc + federal_pct + refundable + max_bft_3 |
                          year + fipsst + birth_year | 0 | cluster,
                        data    = yrbs_all_model,
                        weights = yrbs_all_model$weight_2)

# Get values from models
clust_df <- NULL

# Adolescents (FE only, nested clusters)
clust_df <- make_coef_df(clust_df, twfe_min_sad_1c, "Adolescents (FE only, nested clusters)")
clust_df <- make_coef_df(clust_df, twfe_min_con_1c, "Adolescents (FE only, nested clusters)")
clust_df <- make_coef_df(clust_df, twfe_min_att_1c, "Adolescents (FE only, nested clusters)")
clust_df <- make_coef_df(clust_df, twfe_min_alc_1c, "Adolescents (FE only, nested clusters)")
clust_df <- make_coef_df(clust_df, twfe_min_mjn_1c, "Adolescents (FE only, nested clusters)")
clust_df <- make_coef_df(clust_df, twfe_min_fgh_1c, "Adolescents (FE only, nested clusters)")

# Adolescents (fully adj., nested clusters)
clust_df <- make_coef_df(clust_df, twfe_min_sad_2c, "Adolescents (fully adj., nested clusters)")
clust_df <- make_coef_df(clust_df, twfe_min_con_2c, "Adolescents (fully adj., nested clusters)")
clust_df <- make_coef_df(clust_df, twfe_min_att_2c, "Adolescents (fully adj., nested clusters)")
clust_df <- make_coef_df(clust_df, twfe_min_alc_2c, "Adolescents (fully adj., nested clusters)")
clust_df <- make_coef_df(clust_df, twfe_min_mjn_2c, "Adolescents (fully adj., nested clusters)")
clust_df <- make_coef_df(clust_df, twfe_min_fgh_2c, "Adolescents (fully adj., nested clusters)")

# Add standard TWFE models for comparison
# Adolescents (FE only, state clusters)
clust_df <- make_coef_df(clust_df, twfe_min_sad_1, "Adolescents (FE only, state clusters)")
clust_df <- make_coef_df(clust_df, twfe_min_con_1, "Adolescents (FE only, state clusters)")
clust_df <- make_coef_df(clust_df, twfe_min_att_1, "Adolescents (FE only, state clusters)")
clust_df <- make_coef_df(clust_df, twfe_min_alc_1, "Adolescents (FE only, state clusters)")
clust_df <- make_coef_df(clust_df, twfe_min_mjn_1, "Adolescents (FE only, state clusters)")
clust_df <- make_coef_df(clust_df, twfe_min_fgh_1, "Adolescents (FE only, state clusters)")

# Adolescents (fully adj., state clusters)
clust_df <- make_coef_df(clust_df, twfe_min_sad_2, "Adolescents (fully adj., state clusters)")
clust_df <- make_coef_df(clust_df, twfe_min_con_2, "Adolescents (fully adj., state clusters)")
clust_df <- make_coef_df(clust_df, twfe_min_att_2, "Adolescents (fully adj., state clusters)")
clust_df <- make_coef_df(clust_df, twfe_min_alc_2, "Adolescents (fully adj., state clusters)")
clust_df <- make_coef_df(clust_df, twfe_min_mjn_2, "Adolescents (fully adj., state clusters)")
clust_df <- make_coef_df(clust_df, twfe_min_fgh_2, "Adolescents (fully adj., state clusters)")

# Clean dataframe of coefficients
clust_df <- clean_coef_df(clust_df)

# Get min. and max. N
min(clust_df$n); max(clust_df$n)

# Generate coefficient plot: Nested clusters
plot_clust <- print_coef_plot(
  clust_df,
  Y_TITLE    = "\nAssociation of $1 increase in minimum wage\nwith rates of mental health outcomes",
  Y_MIN      = -0.045,
  Y_MAX      =  0.045
  )

# Adjust legend
plot_clust <- plot_clust +
  guides(shape = guide_legend(reverse=T, nrow=4),
         color = guide_legend(reverse=T, nrow=4))

# Export figure
ggsave(plot=plot_clust, file="Exhibits/YRBS coefficient plot, nested clusters.pdf",
       width=5.75, height=5, units='in', dpi=600)

##############################################################################
# Lifetime minimum wage models
##############################################################################

# Sad or hopeless
life_min_sad_1a <- felm(sad_hopeless ~ wage_life_nom |
                          year + fipsst + birth_year | 0 | fipsst,
                        data    = yrbs_all_model,
                        weights = yrbs_all_model$weight_2)
life_min_sad_2a <- felm(sad_hopeless ~ wage_life_nom +
                          age_2 + sex_2 + race_eth_2 + grade_2 +
                          elig_1_5 + elig_6_18 + has_eitc + federal_pct + refundable + max_bft_3 |
                          year + fipsst + birth_year | 0 | fipsst,
                        data    = yrbs_all_model,
                        weights = yrbs_all_model$weight_2)

life_min_sad_1b <- felm(sad_hopeless ~ wage_life_inf |
                          year + fipsst + birth_year | 0 | fipsst,
                        data    = yrbs_all_model,
                        weights = yrbs_all_model$weight_2)
life_min_sad_2b <- felm(sad_hopeless ~ wage_life_inf +
                          age_2 + sex_2 + race_eth_2 + grade_2 +
                          elig_1_5 + elig_6_18 + has_eitc + federal_pct + refundable + max_bft_3 |
                          year + fipsst + birth_year | 0 | fipsst,
                        data    = yrbs_all_model,
                        weights = yrbs_all_model$weight_2)

# Considered suicide
life_min_con_1a <- felm(cons_suicide ~ wage_life_nom |
                          year + fipsst + birth_year | 0 | fipsst,
                        data    = yrbs_all_model,
                        weights = yrbs_all_model$weight_2)
life_min_con_2a <- felm(cons_suicide ~ wage_life_nom +
                          age_2 + sex_2 + race_eth_2 + grade_2 +
                          elig_1_5 + elig_6_18 + has_eitc + federal_pct + refundable + max_bft_3 |
                          year + fipsst + birth_year | 0 | fipsst,
                        data    = yrbs_all_model,
                        weights = yrbs_all_model$weight_2)

life_min_con_1b <- felm(cons_suicide ~ wage_life_inf |
                          year + fipsst + birth_year | 0 | fipsst,
                        data    = yrbs_all_model,
                        weights = yrbs_all_model$weight_2)
life_min_con_2b <- felm(cons_suicide ~ wage_life_inf +
                          age_2 + sex_2 + race_eth_2 + grade_2 +
                          elig_1_5 + elig_6_18 + has_eitc + federal_pct + refundable + max_bft_3 |
                          year + fipsst + birth_year | 0 | fipsst,
                        data    = yrbs_all_model,
                        weights = yrbs_all_model$weight_2)

# Attempted suicide
life_min_att_1a <- felm(suicide_att ~ wage_life_nom |
                          year + fipsst + birth_year | 0 | fipsst,
                        data    = yrbs_all_model,
                        weights = yrbs_all_model$weight_2)
life_min_att_2a <- felm(suicide_att ~ wage_life_nom +
                          age_2 + sex_2 + race_eth_2 + grade_2 +
                          elig_1_5 + elig_6_18 + has_eitc + federal_pct + refundable + max_bft_3 |
                          year + fipsst + birth_year | 0 | fipsst,
                        data    = yrbs_all_model,
                        weights = yrbs_all_model$weight_2)

life_min_att_1b <- felm(suicide_att ~ wage_life_inf |
                          year + fipsst + birth_year | 0 | fipsst,
                        data    = yrbs_all_model,
                        weights = yrbs_all_model$weight_2)
life_min_att_2b <- felm(suicide_att ~ wage_life_inf +
                          age_2 + sex_2 + race_eth_2 + grade_2 +
                          elig_1_5 + elig_6_18 + has_eitc + federal_pct + refundable + max_bft_3 |
                          year + fipsst + birth_year | 0 | fipsst,
                        data    = yrbs_all_model,
                        weights = yrbs_all_model$weight_2)

# Alcohol use
life_min_alc_1a <- felm(alcohol ~ wage_life_nom |
                          year + fipsst + birth_year | 0 | fipsst,
                        data    = yrbs_all_model,
                        weights = yrbs_all_model$weight_2)
life_min_alc_2a <- felm(alcohol ~ wage_life_nom +
                          age_2 + sex_2 + race_eth_2 + grade_2 +
                          elig_1_5 + elig_6_18 + has_eitc + federal_pct + refundable + max_bft_3 |
                          year + fipsst + birth_year | 0 | fipsst,
                        data    = yrbs_all_model,
                        weights = yrbs_all_model$weight_2)

life_min_alc_1b <- felm(alcohol ~ wage_life_inf |
                          year + fipsst + birth_year | 0 | fipsst,
                        data    = yrbs_all_model,
                        weights = yrbs_all_model$weight_2)
life_min_alc_2b <- felm(alcohol ~ wage_life_inf +
                          age_2 + sex_2 + race_eth_2 + grade_2 +
                          elig_1_5 + elig_6_18 + has_eitc + federal_pct + refundable + max_bft_3 |
                          year + fipsst + birth_year | 0 | fipsst,
                        data    = yrbs_all_model,
                        weights = yrbs_all_model$weight_2)

# Marijuana use
life_min_mjn_1a <- felm(marijuana ~ wage_life_nom |
                          year + fipsst + birth_year | 0 | fipsst,
                        data    = yrbs_all_model,
                        weights = yrbs_all_model$weight_2)
life_min_mjn_2a <- felm(marijuana ~ wage_life_nom +
                          age_2 + sex_2 + race_eth_2 + grade_2 +
                          elig_1_5 + elig_6_18 + has_eitc + federal_pct + refundable + max_bft_3 |
                          year + fipsst + birth_year | 0 | fipsst,
                        data    = yrbs_all_model,
                        weights = yrbs_all_model$weight_2)

life_min_mjn_1b <- felm(marijuana ~ wage_life_inf |
                          year + fipsst + birth_year | 0 | fipsst,
                        data    = yrbs_all_model,
                        weights = yrbs_all_model$weight_2)
life_min_mjn_2b <- felm(marijuana ~ wage_life_inf +
                          age_2 + sex_2 + race_eth_2 + grade_2 +
                          elig_1_5 + elig_6_18 + has_eitc + federal_pct + refundable + max_bft_3 |
                          year + fipsst + birth_year | 0 | fipsst,
                        data    = yrbs_all_model,
                        weights = yrbs_all_model$weight_2)

# Physical fight
life_min_fgh_1a <- felm(fight ~ wage_life_nom |
                          year + fipsst + birth_year | 0 | fipsst,
                        data    = yrbs_all_model,
                        weights = yrbs_all_model$weight_2)
life_min_fgh_2a <- felm(fight ~ wage_life_nom +
                          age_2 + sex_2 + race_eth_2 + grade_2 +
                          elig_1_5 + elig_6_18 + has_eitc + federal_pct + refundable + max_bft_3 |
                          year + fipsst + birth_year | 0 | fipsst,
                        data    = yrbs_all_model,
                        weights = yrbs_all_model$weight_2)

life_min_fgh_1b <- felm(fight ~ wage_life_inf |
                          year + fipsst + birth_year | 0 | fipsst,
                        data    = yrbs_all_model,
                        weights = yrbs_all_model$weight_2)
life_min_fgh_2b <- felm(fight ~ wage_life_inf +
                          age_2 + sex_2 + race_eth_2 + grade_2 +
                          elig_1_5 + elig_6_18 + has_eitc + federal_pct + refundable + max_bft_3 |
                          year + fipsst + birth_year | 0 | fipsst,
                        data    = yrbs_all_model,
                        weights = yrbs_all_model$weight_2)

# Get values from models
life_df <- NULL

# Adolescents (FE only)
life_df <- make_coef_df(life_df, life_min_sad_1a, "Adolescents (FE only), nominal wages")
life_df <- make_coef_df(life_df, life_min_con_1a, "Adolescents (FE only), nominal wages")
life_df <- make_coef_df(life_df, life_min_att_1a, "Adolescents (FE only), nominal wages")
life_df <- make_coef_df(life_df, life_min_alc_1a, "Adolescents (FE only), nominal wages")
life_df <- make_coef_df(life_df, life_min_mjn_1a, "Adolescents (FE only), nominal wages")
life_df <- make_coef_df(life_df, life_min_fgh_1a, "Adolescents (FE only), nominal wages")

life_df <- make_coef_df(life_df, life_min_sad_1b, "Adolescents (FE only), real wages")
life_df <- make_coef_df(life_df, life_min_con_1b, "Adolescents (FE only), real wages")
life_df <- make_coef_df(life_df, life_min_att_1b, "Adolescents (FE only), real wages")
life_df <- make_coef_df(life_df, life_min_alc_1b, "Adolescents (FE only), real wages")
life_df <- make_coef_df(life_df, life_min_mjn_1b, "Adolescents (FE only), real wages")
life_df <- make_coef_df(life_df, life_min_fgh_1b, "Adolescents (FE only), real wages")

# Adolescents (fully adj.)
life_df <- make_coef_df(life_df, life_min_sad_2a, "Adolescents (fully adj.), nominal wages")
life_df <- make_coef_df(life_df, life_min_con_2a, "Adolescents (fully adj.), nominal wages")
life_df <- make_coef_df(life_df, life_min_att_2a, "Adolescents (fully adj.), nominal wages")
life_df <- make_coef_df(life_df, life_min_alc_2a, "Adolescents (fully adj.), nominal wages")
life_df <- make_coef_df(life_df, life_min_mjn_2a, "Adolescents (fully adj.), nominal wages")
life_df <- make_coef_df(life_df, life_min_fgh_2a, "Adolescents (fully adj.), nominal wages")

life_df <- make_coef_df(life_df, life_min_sad_2b, "Adolescents (fully adj.), real wages")
life_df <- make_coef_df(life_df, life_min_con_2b, "Adolescents (fully adj.), real wages")
life_df <- make_coef_df(life_df, life_min_att_2b, "Adolescents (fully adj.), real wages")
life_df <- make_coef_df(life_df, life_min_alc_2b, "Adolescents (fully adj.), real wages")
life_df <- make_coef_df(life_df, life_min_mjn_2b, "Adolescents (fully adj.), real wages")
life_df <- make_coef_df(life_df, life_min_fgh_2b, "Adolescents (fully adj.), real wages")

# Clean dataframe of coefficients
life_df <- clean_coef_df(life_df)

# Get min. and max. N
min(life_df$n); max(life_df$n)

# Generate coefficient plot: Lifetime wage
plot_life <- print_coef_plot(
  life_df,
  Y_TITLE    = "\nAssociation of $1 increase in lifetime minimum\nwage with rates of mental health outcomes",
  Y_MIN      = -0.045,
  Y_MAX      =  0.045
  )

# Adjust legend
plot_life <- plot_life +
  guides(shape = guide_legend(reverse=T, nrow=4),
         color = guide_legend(reverse=T, nrow=4))

# Export figure
ggsave(plot=plot_life, file="Exhibits/YRBS coefficient plot, lifetime.pdf",
       width=5.75, height=5, units='in', dpi=600)

##############################################################################
# Preparation for DID/event study models
##############################################################################

# Treat year as categorical
yrbs_all_model$year_cat <- relevel(as.factor(yrbs_all_model$year), "2013")

# Define treated and control states
TREATED <- c("AR","DE","HI","MD","MT","NE","NJ","NY","SD","WV")
CONTROL <- c("AL","GA","IA","ID","IN","KS","KY","LA","MS","NC",
             "ND","NH","OK","PA","SC","TN","TX","UT","WI","WY")
TREATED_2 <- cdlTools::fips(TREATED, to="FIPS")
CONTROL_2 <- cdlTools::fips(CONTROL, to="FIPS")

# Code treatment groups
yrbs_all_model <- yrbs_all_model %>% mutate(
  event_treated = case_when(
    fipsst %in% TREATED_2 ~ 1,
    fipsst %in% CONTROL_2 ~ 0
  ))

# Code treatment years
yrbs_all_model <- yrbs_all_model %>% mutate(
  treated_years = case_when(
    year %in% c(2011:2013) ~ 0,
    year %in% c(2014:2021) ~ 1
  ))

# Subset dataset to relevant years
yrbs_event <- subset(yrbs_all_model, year %in% c(2011:2021) &
                       fipsst %in% c(TREATED_2, CONTROL_2))

# Function to extract and clean values from event study models.
# Requires: Df for saving coefficients, "lfe" model, title of model.
# Returns: Df of values, requires cleaning before passing to "print_event_plot".
make_event_df <- function(event_df, event_model, TITLE) {
  
  # Get outcome labels
  if (event_model$lhs == "sad_hopeless") {
    outcome  <- "Sad or hopeless in past year"}
  
  if (event_model$lhs == "cons_suicide") {
    outcome  <- "Considered suicide in past year"}
  
  if (event_model$lhs == "suicide_att") {
    outcome  <- "Attempted suicide in past year"}
  
  if (event_model$lhs == "alcohol") {
    outcome  <- "Alcohol use in past month"}
  
  if (event_model$lhs == "marijuana") {
    outcome  <- "Marijuana use in past month"}
  
  if (event_model$lhs == "fight") {
    outcome  <- "Physical fight in past year"}
  
  # Get coefficients of model
  coef <- data.frame(cbind(outcome, event_model$coefficients, event_model$cse, TITLE))
  colnames(coef) <- c("outcome", "effect", "se", "name"); coef$var <- rownames(coef)
  
  # Only include year coefficients
  coef <- coef %>% filter(grepl("event_treated:year_cat", var))
  
  # Rename period variable
  coef$period <- gsub("event_treated:year_cat", "", coef$var)
  
  # Add zero period
  coef <- data.frame(rbind(coef, c(
    outcome, 0, 0, TITLE, "event_treated:year_cat2013", 2013)))
  
  # Treat columns as numeric
  coef$effect  <- as.numeric(coef$effect)
  coef$se      <- as.numeric(coef$se)
  coef$period  <- as.numeric(coef$period)
  
  # Add to existing event study df
  event_df <- rbind(event_df, coef)
  
  return(event_df)
}

# Function to generate event study plots.
# Requires: Df of coefficients, y title, and y limits.
# Returns: Formatted event study plot.
print_event_plot <- function(event_df, Y_TITLE, Y_MIN, Y_MAX) {
  
  # Generate event study plot
  plot <- ggplot(event_df, aes(x=period, y=effect, group=name, color=name)) +
    
    # Null line
    geom_hline(yintercept=0, color="black") +
    
    # Event line
    geom_vline(xintercept=2013.5, color="black", linetype="dashed") +
    
    # Confidence intervals
    geom_errorbar(aes(ymin=effect-1.96*se, ymax=effect+1.96*se),
                  position = position_dodge(width=0.6), width=0, linewidth=0.8) +
    geom_errorbar(aes(ymin=effect-2.64*se, ymax=effect+2.64*se),
                  position = position_dodge(width=0.6), width=0.3, alpha=0.5) +
    
    # Lines and point estimates with distinct shapes
    geom_line(position = position_dodge(width = 0.5)) +
    geom_point(position = position_dodge(width = 0.5), size=1, aes(shape=name)) +
    scale_shape_manual(values = 1:nlevels(event_df$name)) +
    
    # Titles
    xlab("Year") +
    ylab(Y_TITLE) +
    
    # Theme modifications
    theme_test() +
    theme(legend.position = "bottom",
          text = element_text(size = 10, face = "bold"),
          axis.title.x = element_blank(),
          axis.ticks = element_blank(),
          strip.background = element_blank(),
          legend.title = element_blank(),
          panel.grid.major.x = element_blank(),
          panel.grid.major.y = element_line(color="light gray", linewidth=0.5),
          panel.grid.minor.y = element_line(color="light gray", linewidth=0.25)) +
    scale_y_continuous(breaks=seq(-0.1, 0.1, 0.05),
                       minor_breaks = seq(-0.1, 0.1, 0.025),
                       limits=c(Y_MIN, Y_MAX),
                       labels = function(x) paste0(x*100," pp")) +
    scale_x_continuous(breaks=seq(2011,2021,2)) +
    annotate("text", x=2016.5, y=0.075, label="Raised minimum wage",
             fontface=2, size=2.5, hjust=0.5, vjust=0.5) +
    annotate("rect", xmin=2013.5, xmax=2021.5,
             ymin=-Inf, ymax=Inf, alpha=0.2, fill="grey70") +
    scale_color_grey(start=0.6, end=0) +
    facet_wrap(.~outcome, ncol=2)
  
  # Return plot
  return(plot)
}

##############################################################################
# Graphs: Descriptives for DID/event studies
##############################################################################

# Make dataframe of states
library(datasets)
event_map_df <- data.frame(datasets::state.abb)
colnames(event_map_df) <- "state_abb"
event_map_df <- rbind(event_map_df, "DC")

# Add FIPS codes
event_map_df$fips <- cdlTools::fips(event_map_df$state_abb, to="FIPS")

# Code treatment groups
event_map_df <- event_map_df %>% mutate(
  value = case_when(
    state_abb %in% TREATED ~ "Raised min. wage in 2014-2015",
    state_abb %in% CONTROL ~ "Used federal min. from 2011-2021",
    TRUE ~ "Not included"
  ))
event_map_df$value <- factor(event_map_df$value, levels = c("Raised min. wage in 2014-2015", "Used federal min. from 2011-2021", "Not included"))

# Map for event study states
event_map <- plot_usmap(regions="states", data=event_map_df,
                        values="value", size=0.4) +
  theme(legend.position = "right",
        text = element_text(size = 10, face = "bold"),
        plot.title = element_text(size=14, face="bold", hjust=0.5)) +
  scale_fill_manual(name = "", values = c("grey30","grey60","grey99"))

# Subset minimum wage df to treated state-years
min_wage_event <- subset(min_wage_long, fipsst %in% c(CONTROL_2, TREATED_2) &
                           year %in% c(2011:2021))

# Graph minimum wages of treated states
event_desc_plot <- ggplot(subset(min_wage_event, fipsst %in% TREATED_2),
                          aes(x=year, y=effective_min_wage_l, group=State.or.other.jurisdiction)) +
  geom_hline(yintercept=7.25, color="black") +
  geom_vline(xintercept=2013, color="black", linetype="dashed") +
  geom_vline(xintercept=2015, color="black", linetype="dashed") +
  geom_line(aes(group=State.or.other.jurisdiction)) +
  ylab("Effective minimum wage") +
  xlab("Year") +
  theme_test() +
  theme(legend.position = "bottom",
        text = element_text(size = 10, face = "bold"),
        axis.title.x = element_blank(),
        axis.text.x = element_text(angle=45, hjust=1, vjust=1),
        axis.ticks = element_blank(),
        strip.background = element_blank(),
        legend.title = element_blank(),
        panel.grid.major.x = element_blank(),
        panel.grid.major.y = element_line(color="light gray", linewidth=0.5),
        panel.grid.minor.y = element_line(color="light gray", linewidth=0.25)) +
  scale_y_continuous(limits = c(6,12.5),
                     breaks = seq(0, 15, 2),
                     minor_breaks = seq(0, 15, 0.5),
                     labels = function(x) paste0("$",x)) +
  scale_x_continuous(limits = c(2011, 2021),
                     breaks = seq(2011, 2021, 2)) +
  facet_wrap(.~State.or.other.jurisdiction, nrow=2)

# Compile figures
event_desc_1 <- plot_grid(event_map, event_desc_plot, rel_heights=c(0.7,1),
                          nrow=2, labels=c("A.", "B."))

# Export figure
ggsave(plot=event_desc_1, file="Exhibits/YRBS event studies, descriptives 1.pdf",
       width=6, height=5, units='in', dpi=600)

# Get mean wages in treated states
yrbs_all_model %>%
  filter(fipsst %in% TREATED_2) %>%
  group_by(treated_years) %>%
  summarise(mean = weighted.mean(x=effective_min_wage_l, w=weight_2))

# Get mean wages in treated states by year
mean_tx_wages <- yrbs_all_model %>%
  filter(fipsst %in% TREATED_2) %>%
  group_by(year) %>%
  summarise(mean = weighted.mean(x=effective_min_wage_l, w=weight_2))
mean_tx_wages$year <- as.numeric(as.character(mean_tx_wages$year))
mean_tx_wages$mean <- as.numeric(as.character(mean_tx_wages$mean))

# Graph minimum wages of treated states
mean_wage_plot <- ggplot(mean_tx_wages, aes(x=year, y=mean)) +
  geom_hline(yintercept=7.25, color="black") +
  geom_vline(xintercept=2013, color="black", linetype="dashed") +
  geom_vline(xintercept=2015, color="black", linetype="dashed") +
  geom_point() + geom_line() +
  ylab("Weighted mean of effective\nmin. wage in treated states") +
  xlab("Year") +
  theme_test() +
  theme(legend.position = "bottom",
        text = element_text(size = 10, face = "bold"),
        axis.title.x = element_blank(),
        axis.ticks = element_blank(),
        strip.background = element_blank(),
        panel.grid.major.x = element_blank(),
        panel.grid.major.y = element_line(color="light gray", linewidth=0.5),
        panel.grid.minor.y = element_line(color="light gray", linewidth=0.25)) +
  scale_y_continuous(limits = c(6,12.5),
                     breaks = seq(0, 15, 2),
                     minor_breaks = seq(0, 15, 0.5),
                     labels = function(x) paste0("$",x)) +
  scale_x_continuous(limits = c(2010.5, 2021),
                     breaks = seq(2011, 2021, 2)) +
  geom_text(aes(label = paste0("$", format(round(mean, 2), nsmall = 2))),
            nudge_x=-0.3, nudge_y=0.4, size=3, fontface="bold")

# Export figure
ggsave(plot=mean_wage_plot, file="Exhibits/YRBS event studies, descriptives 2.pdf",
       width=4, height=3, units='in', dpi=600)

##############################################################################
# Main difference-in-differences models
##############################################################################

# Sad or hopeless
did_sad_1 <- felm(sad_hopeless ~ event_treated*treated_years |
                    year + fipsst + birth_year | 0 | fipsst,
                  data    = yrbs_event,
                  weights = yrbs_event$weight_2)
did_sad_2 <- felm(sad_hopeless ~ event_treated*treated_years +
                    age_2 + sex_2 + race_eth_2 + grade_2 +
                    elig_1_5 + elig_6_18 + has_eitc + federal_pct + refundable + max_bft_3 |
                    year + fipsst + birth_year | 0 | fipsst,
                  data    = yrbs_event,
                  weights = yrbs_event$weight_2)

# Considered suicide
did_con_1 <- felm(cons_suicide ~ event_treated*treated_years |
                    year + fipsst + birth_year | 0 | fipsst,
                  data    = yrbs_event,
                  weights = yrbs_event$weight_2)
did_con_2 <- felm(cons_suicide ~ event_treated*treated_years +
                    age_2 + sex_2 + race_eth_2 + grade_2 +
                    elig_1_5 + elig_6_18 + has_eitc + federal_pct + refundable + max_bft_3 |
                    year + fipsst + birth_year | 0 | fipsst,
                  data    = yrbs_event,
                  weights = yrbs_event$weight_2)

# Attempted suicide
did_att_1 <- felm(suicide_att ~ event_treated*treated_years |
                    year + fipsst + birth_year | 0 | fipsst,
                  data    = yrbs_event,
                  weights = yrbs_event$weight_2)
did_att_2 <- felm(suicide_att ~ event_treated*treated_years +
                    age_2 + sex_2 + race_eth_2 + grade_2 +
                    elig_1_5 + elig_6_18 + has_eitc + federal_pct + refundable + max_bft_3 |
                    year + fipsst + birth_year | 0 | fipsst,
                  data    = yrbs_event,
                  weights = yrbs_event$weight_2)

# Alcohol use
did_alc_1 <- felm(alcohol ~ event_treated*treated_years |
                    year + fipsst + birth_year | 0 | fipsst,
                  data    = yrbs_event,
                  weights = yrbs_event$weight_2)
did_alc_2 <- felm(alcohol ~ event_treated*treated_years +
                    age_2 + sex_2 + race_eth_2 + grade_2 +
                    elig_1_5 + elig_6_18 + has_eitc + federal_pct + refundable + max_bft_3 |
                    year + fipsst + birth_year | 0 | fipsst,
                  data    = yrbs_event,
                  weights = yrbs_event$weight_2)

# Alcohol use
did_mjn_1 <- felm(marijuana ~ event_treated*treated_years |
                    year + fipsst + birth_year | 0 | fipsst,
                  data    = yrbs_event,
                  weights = yrbs_event$weight_2)
did_mjn_2 <- felm(marijuana ~ event_treated*treated_years +
                    age_2 + sex_2 + race_eth_2 + grade_2 +
                    elig_1_5 + elig_6_18 + has_eitc + federal_pct + refundable + max_bft_3 |
                    year + fipsst + birth_year | 0 | fipsst,
                  data    = yrbs_event,
                  weights = yrbs_event$weight_2)

# Physical fight
did_fgh_1 <- felm(fight ~ event_treated*treated_years |
                    year + fipsst + birth_year | 0 | fipsst,
                  data    = yrbs_event,
                  weights = yrbs_event$weight_2)
did_fgh_2 <- felm(fight ~ event_treated*treated_years +
                    age_2 + sex_2 + race_eth_2 + grade_2 +
                    elig_1_5 + elig_6_18 + has_eitc + federal_pct + refundable + max_bft_3 |
                    year + fipsst + birth_year | 0 | fipsst,
                  data    = yrbs_event,
                  weights = yrbs_event$weight_2)

# Additional row in table
cov_row <- as.data.frame(
  rbind(cbind("Demographic controls",     "No",  "Yes", "No",  "Yes", "No",  "Yes"),
        cbind("State policy controls",    "No",  "Yes", "No",  "Yes", "No",  "Yes"),
        cbind("State & age-by-year FEs",  "Yes", "Yes", "Yes", "Yes", "Yes", "Yes"),
        cbind("Cluster robust SEs", "State", "State", "State", "State", "State", "State")))

# Compile results into tables
modelsummary(list("(1)" = did_sad_1,
                  "(2)" = did_sad_2,
                  "(1)" = did_con_1,
                  "(2)" = did_con_2,
                  "(1)" = did_att_1,
                  "(2)" = did_att_2),
             gof_omit    = "Log*|AIC|BIC|F|RMSE|Std",
             coef_omit   = "^(?!event_treated:treated_years)",
             coef_rename = c("event_treated:treated_years" =
                               "Effect of raise in wage"),
             statistic   = c("[{conf.low}, {conf.high}]"),
             # conf_level  = 0.99167,
             add_rows    = cov_row,
             fmt         = fmt_statistic("estimate"=4, "conf.low"=4, "conf.high"=4)) %>%
  add_header_above(c(" " = 1, "Sad or hopeless" = 2, "Cons. suicide" = 2, "Att. suicide" = 2))

modelsummary(list("(1)" = did_alc_1,
                  "(2)" = did_alc_2,
                  "(1)" = did_mjn_1,
                  "(2)" = did_mjn_2,
                  "(1)" = did_fgh_1,
                  "(2)" = did_fgh_2),
             gof_omit    = "Log*|AIC|BIC|F|RMSE|Std",
             coef_omit   = "^(?!event_treated:treated_years)",
             coef_rename = c("event_treated:treated_years" =
                               "Effect of raise in wage"),
             statistic   = c("[{conf.low}, {conf.high}]"),
             # conf_level  = 0.99167,
             add_rows    = cov_row,
             fmt         = fmt_statistic("estimate"=4, "conf.low"=4, "conf.high"=4)) %>%
  add_header_above(c(" " = 1, "Alcohol" = 2, "Marijuana" = 2, "Phys. fight" = 2))

##############################################################################
# Main event studies
##############################################################################

# Sad or hopeless
event_sad_1 <- felm(sad_hopeless ~ event_treated*year_cat |
                      year + fipsst + birth_year | 0 | fipsst,
                    data    = yrbs_event,
                    weights = yrbs_event$weight_2)
event_sad_2 <- felm(sad_hopeless ~ event_treated*year_cat +
                      age_2 + sex_2 + race_eth_2 + grade_2 +
                      elig_1_5 + elig_6_18 + has_eitc + federal_pct + refundable + max_bft_3 |
                      year + fipsst + birth_year | 0 | fipsst,
                    data    = yrbs_event,
                    weights = yrbs_event$weight_2)

# Considered suicide
event_con_1 <- felm(cons_suicide ~ event_treated*year_cat |
                      year + fipsst + birth_year | 0 | fipsst,
                    data    = yrbs_event,
                    weights = yrbs_event$weight_2)
event_con_2 <- felm(cons_suicide ~ event_treated*year_cat +
                      age_2 + sex_2 + race_eth_2 + grade_2 +
                      elig_1_5 + elig_6_18 + has_eitc + federal_pct + refundable + max_bft_3 |
                      year + fipsst + birth_year | 0 | fipsst,
                    data    = yrbs_event,
                    weights = yrbs_event$weight_2)

# Attempted suicide
event_att_1 <- felm(suicide_att ~ event_treated*year_cat |
                      year + fipsst + birth_year | 0 | fipsst,
                    data    = yrbs_event,
                    weights = yrbs_event$weight_2)
event_att_2 <- felm(suicide_att ~ event_treated*year_cat +
                      age_2 + sex_2 + race_eth_2 + grade_2 +
                      elig_1_5 + elig_6_18 + has_eitc + federal_pct + refundable + max_bft_3 |
                      year + fipsst + birth_year | 0 | fipsst,
                    data    = yrbs_event,
                    weights = yrbs_event$weight_2)

# Alcohol use
event_alc_1 <- felm(alcohol ~ event_treated*year_cat |
                      year + fipsst + birth_year | 0 | fipsst,
                    data    = yrbs_event,
                    weights = yrbs_event$weight_2)
event_alc_2 <- felm(alcohol ~ event_treated*year_cat +
                      age_2 + sex_2 + race_eth_2 + grade_2 +
                      elig_1_5 + elig_6_18 + has_eitc + federal_pct + refundable + max_bft_3 |
                      year + fipsst + birth_year | 0 | fipsst,
                    data    = yrbs_event,
                    weights = yrbs_event$weight_2)

# Alcohol use
event_mjn_1 <- felm(marijuana ~ event_treated*year_cat |
                      year + fipsst + birth_year | 0 | fipsst,
                    data    = yrbs_event,
                    weights = yrbs_event$weight_2)
event_mjn_2 <- felm(marijuana ~ event_treated*year_cat +
                      age_2 + sex_2 + race_eth_2 + grade_2 +
                      elig_1_5 + elig_6_18 + has_eitc + federal_pct + refundable + max_bft_3 |
                      year + fipsst + birth_year | 0 | fipsst,
                    data    = yrbs_event,
                    weights = yrbs_event$weight_2)

# Physical fight
event_fgh_1 <- felm(fight ~ event_treated*year_cat |
                      year + fipsst + birth_year | 0 | fipsst,
                    data    = yrbs_event,
                    weights = yrbs_event$weight_2)
event_fgh_2 <- felm(fight ~ event_treated*year_cat +
                      age_2 + sex_2 + race_eth_2 + grade_2 +
                      elig_1_5 + elig_6_18 + has_eitc + federal_pct + refundable + max_bft_3 |
                      year + fipsst + birth_year | 0 | fipsst,
                    data    = yrbs_event,
                    weights = yrbs_event$weight_2)

# Get values from models
event_df <- NULL

# Adolescents (FE only)
event_df <- make_event_df(event_df, event_sad_1, "Adolescents (FE only)")
event_df <- make_event_df(event_df, event_con_1, "Adolescents (FE only)")
event_df <- make_event_df(event_df, event_att_1, "Adolescents (FE only)")
event_df <- make_event_df(event_df, event_alc_1, "Adolescents (FE only)")
event_df <- make_event_df(event_df, event_mjn_1, "Adolescents (FE only)")
event_df <- make_event_df(event_df, event_fgh_1, "Adolescents (FE only)")

# Adolescents (fully adjusted)
event_df <- make_event_df(event_df, event_sad_2, "Adolescents (fully adjusted)")
event_df <- make_event_df(event_df, event_con_2, "Adolescents (fully adjusted)")
event_df <- make_event_df(event_df, event_att_2, "Adolescents (fully adjusted)")
event_df <- make_event_df(event_df, event_alc_2, "Adolescents (fully adjusted)")
event_df <- make_event_df(event_df, event_mjn_2, "Adolescents (fully adjusted)")
event_df <- make_event_df(event_df, event_fgh_2, "Adolescents (fully adjusted)")

# Reorder outcomes
event_df$outcome <- factor(
  event_df$outcome, levels=c("Sad or hopeless in past year", "Considered suicide in past year", "Attempted suicide in past year", "Alcohol use in past month", "Marijuana use in past month", "Physical fight in past year"))

# Get min. and max. N
event_n <-
  (c(length(event_sad_1$residuals), length(event_sad_2$residuals),
     length(event_con_1$residuals), length(event_con_2$residuals),
     length(event_att_1$residuals), length(event_att_2$residuals),
     length(event_alc_1$residuals), length(event_alc_2$residuals),
     length(event_mjn_1$residuals), length(event_mjn_2$residuals),
     length(event_fgh_1$residuals), length(event_fgh_2$residuals)))
min(event_n); max(event_n)

# Generate event study plot: Main
event_plot <- print_event_plot(
  event_df,
  Y_TITLE    = "Effect of raising minimum wage\non indicated outcome by year",
  Y_MIN      = -0.1,
  Y_MAX      =  0.1
)

# Export figure
ggsave(plot=event_plot, file="Exhibits/YRBS event studies, main.pdf",
       width=6, height=6, units='in', dpi=600)

##############################################################################
# Event study robustness check: Nested clusters
##############################################################################

# Sad or hopeless
event_sad_1c <- felm(sad_hopeless ~ event_treated*year_cat |
                       year + fipsst + birth_year | 0 | cluster,
                     data    = yrbs_event,
                     weights = yrbs_event$weight_2)
event_sad_2c <- felm(sad_hopeless ~ event_treated*year_cat +
                       age_2 + sex_2 + race_eth_2 + grade_2 +
                       elig_1_5 + elig_6_18 + has_eitc + federal_pct + refundable + max_bft_3 |
                       year + fipsst + birth_year | 0 | cluster,
                     data    = yrbs_event,
                     weights = yrbs_event$weight_2)

# Considered suicide
event_con_1c <- felm(cons_suicide ~ event_treated*year_cat |
                       year + fipsst + birth_year | 0 | cluster,
                     data    = yrbs_event,
                     weights = yrbs_event$weight_2)
event_con_2c <- felm(cons_suicide ~ event_treated*year_cat +
                       age_2 + sex_2 + race_eth_2 + grade_2 +
                       elig_1_5 + elig_6_18 + has_eitc + federal_pct + refundable + max_bft_3 |
                       year + fipsst + birth_year | 0 | cluster,
                     data    = yrbs_event,
                     weights = yrbs_event$weight_2)

# Attempted suicide
event_att_1c <- felm(suicide_att ~ event_treated*year_cat |
                       year + fipsst + birth_year | 0 | cluster,
                     data    = yrbs_event,
                     weights = yrbs_event$weight_2)
event_att_2c <- felm(suicide_att ~ event_treated*year_cat +
                       age_2 + sex_2 + race_eth_2 + grade_2 +
                       elig_1_5 + elig_6_18 + has_eitc + federal_pct + refundable + max_bft_3 |
                       year + fipsst + birth_year | 0 | cluster,
                     data    = yrbs_event,
                     weights = yrbs_event$weight_2)

# Alcohol use
event_alc_1c <- felm(alcohol ~ event_treated*year_cat |
                       year + fipsst + birth_year | 0 | cluster,
                     data    = yrbs_event,
                     weights = yrbs_event$weight_2)
event_alc_2c <- felm(alcohol ~ event_treated*year_cat +
                       age_2 + sex_2 + race_eth_2 + grade_2 +
                       elig_1_5 + elig_6_18 + has_eitc + federal_pct + refundable + max_bft_3 |
                       year + fipsst + birth_year | 0 | cluster,
                     data    = yrbs_event,
                     weights = yrbs_event$weight_2)

# Alcohol use
event_mjn_1c <- felm(marijuana ~ event_treated*year_cat |
                       year + fipsst + birth_year | 0 | cluster,
                     data    = yrbs_event,
                     weights = yrbs_event$weight_2)
event_mjn_2c <- felm(marijuana ~ event_treated*year_cat +
                       age_2 + sex_2 + race_eth_2 + grade_2 +
                       elig_1_5 + elig_6_18 + has_eitc + federal_pct + refundable + max_bft_3 |
                       year + fipsst + birth_year | 0 | cluster,
                     data    = yrbs_event,
                     weights = yrbs_event$weight_2)

# Physical fight
event_fgh_1c <- felm(fight ~ event_treated*year_cat |
                       year + fipsst + birth_year | 0 | cluster,
                     data    = yrbs_event,
                     weights = yrbs_event$weight_2)
event_fgh_2c <- felm(fight ~ event_treated*year_cat +
                       age_2 + sex_2 + race_eth_2 + grade_2 +
                       elig_1_5 + elig_6_18 + has_eitc + federal_pct + refundable + max_bft_3 |
                       year + fipsst + birth_year | 0 | cluster,
                     data    = yrbs_event,
                     weights = yrbs_event$weight_2)

# Get values from models
ev_clust_df <- NULL

# Adolescents (FE only, nested clusters)
ev_clust_df <- make_event_df(ev_clust_df, event_sad_1c, "Adolescents (FE only, nested clusters)")
ev_clust_df <- make_event_df(ev_clust_df, event_con_1c, "Adolescents (FE only, nested clusters)")
ev_clust_df <- make_event_df(ev_clust_df, event_att_1c, "Adolescents (FE only, nested clusters)")
ev_clust_df <- make_event_df(ev_clust_df, event_alc_1c, "Adolescents (FE only, nested clusters)")
ev_clust_df <- make_event_df(ev_clust_df, event_mjn_1c, "Adolescents (FE only, nested clusters)")
ev_clust_df <- make_event_df(ev_clust_df, event_fgh_1c, "Adolescents (FE only, nested clusters)")

# Adolescents (fully adj., nested clusters)
ev_clust_df <- make_event_df(ev_clust_df, event_sad_2c,
                             "Adolescents (fully adj., nested clusters)")
ev_clust_df <- make_event_df(ev_clust_df, event_con_2c,
                             "Adolescents (fully adj., nested clusters)")
ev_clust_df <- make_event_df(ev_clust_df, event_att_2c,
                             "Adolescents (fully adj., nested clusters)")
ev_clust_df <- make_event_df(ev_clust_df, event_alc_2c,
                             "Adolescents (fully adj., nested clusters)")
ev_clust_df <- make_event_df(ev_clust_df, event_mjn_2c,
                             "Adolescents (fully adj., nested clusters)")
ev_clust_df <- make_event_df(ev_clust_df, event_fgh_2c,
                             "Adolescents (fully adj., nested clusters)")

# Reorder outcomes
ev_clust_df$outcome <- factor(
  ev_clust_df$outcome, levels=c("Sad or hopeless in past year", "Considered suicide in past year", "Attempted suicide in past year", "Alcohol use in past month", "Marijuana use in past month", "Physical fight in past year"))

# Get min. and max. N
cluster_n <-
  (c(length(event_sad_1c$residuals), length(event_sad_2c$residuals),
     length(event_con_1c$residuals), length(event_con_2c$residuals),
     length(event_att_1c$residuals), length(event_att_2c$residuals),
     length(event_alc_1c$residuals), length(event_alc_2c$residuals),
     length(event_mjn_1c$residuals), length(event_mjn_2c$residuals),
     length(event_fgh_1c$residuals), length(event_fgh_2c$residuals)))
min(cluster_n); max(cluster_n)

# Generate event study plot: Nested clusters
event_clust_plot <- print_event_plot(
  ev_clust_df,
  Y_TITLE    = "Effect of raising minimum wage\non indicated outcome by year",
  Y_MIN      = -0.1,
  Y_MAX      =  0.1
)

# Export figure
ggsave(plot=event_clust_plot, file="Exhibits/YRBS event studies, nested clusters.pdf",
       width=6, height=6, units='in', dpi=600)

##############################################################################
# Event study robustness check: Strictly balanced panels
##############################################################################

# Function to extract and clean values from TWFE models.
# Requires: Df for saving coefficients, "lfe" model, title of model.
# Returns: Df of values, ready to pass to "clean_coef_df".
make_event_df_2 <- function(event_df, event_model, TITLE) {
  
  # Get outcome labels
  if (event_model$lhs == "sad_hopeless") {
    outcome  <- "Sad or hopeless in past year"}
  
  if (event_model$lhs == "cons_suicide") {
    outcome  <- "Considered suicide in past year"}
  
  if (event_model$lhs == "suicide_att") {
    outcome  <- "Attempted suicide in past year"}
  
  if (event_model$lhs == "alcohol") {
    outcome  <- "Alcohol use in past month"}
  
  if (event_model$lhs == "marijuana") {
    outcome  <- "Marijuana use in past month"}
  
  if (event_model$lhs == "fight") {
    outcome  <- "Physical fight in past year"}
  
  # Get coefficients of model
  coef <- data.frame(cbind(outcome, event_model$coefficients, event_model$cse, TITLE))
  colnames(coef) <- c("outcome", "effect", "se", "name"); coef$var <- rownames(coef)
  
  # Only include year coefficients
  coef <- coef %>% filter(grepl("event_tx_balance*.", var))
  
  # Rename period variable
  coef$period <- gsub(".*year_cat", "", coef$var)
  
  # Add zero period
  coef <- data.frame(rbind(coef, c(
    outcome, 0, 0, TITLE, "event_tx_balance:year_cat2013", 2013)))
  
  # Treat columns as numeric
  coef$effect  <- as.numeric(coef$effect)
  coef$se      <- as.numeric(coef$se)
  coef$period  <- as.numeric(coef$period)
  
  # Add to existing event study df
  event_df <- rbind(event_df, coef)
  
  return(event_df)
}

# Define treated and control states
# Note: Balanced panel is different for each outcome
TREATED_B_SAD <- c("AR","HI","MD","MT","NE","NY","WV")
CONTROL_B_SAD <- c("ID","KY","NC","ND","NH","OK","SC","TN")
TREATED_B_CON <- c("AR","HI","MD","MT","NE","NY","WV")
CONTROL_B_CON <- c("ID","KY","NC","ND","NH","OK","TN")
TREATED_B_ATT <- c("AR","HI","MT","NE","NY","WV")
CONTROL_B_ATT <- c("ID","KY","ND","NH","OK","TN")
TREATED_B_ALC <- c("AR","HI","MD","MT","NE","NY","WV")
CONTROL_B_ALC <- c("ID","KY","NC","ND","NH","OK","SC")
TREATED_B_MJN <- c("AR","HI","MD","MT","NE","NY","WV")
CONTROL_B_MJN <- c("ID","KY","NC","ND","NH","OK","SC")
TREATED_B_FGH <- c("AR","HI","MT","NE","NY","WV")
CONTROL_B_FGH <- c("NC","SC")

# Code treatment groups
yrbs_event <- yrbs_event %>% mutate(
  event_tx_balance_sad = case_when(
    fipsst %in% cdlTools::fips(TREATED_B_SAD, to="FIPS") ~ 1,
    fipsst %in% cdlTools::fips(CONTROL_B_SAD, to="FIPS") ~ 0),
  event_tx_balance_con = case_when(
    fipsst %in% cdlTools::fips(TREATED_B_CON, to="FIPS") ~ 1,
    fipsst %in% cdlTools::fips(CONTROL_B_CON, to="FIPS") ~ 0),
  event_tx_balance_att = case_when(
    fipsst %in% cdlTools::fips(TREATED_B_ATT, to="FIPS") ~ 1,
    fipsst %in% cdlTools::fips(CONTROL_B_ATT, to="FIPS") ~ 0),
  event_tx_balance_alc = case_when(
    fipsst %in% cdlTools::fips(TREATED_B_ALC, to="FIPS") ~ 1,
    fipsst %in% cdlTools::fips(CONTROL_B_ALC, to="FIPS") ~ 0),
  event_tx_balance_mjn= case_when(
    fipsst %in% cdlTools::fips(TREATED_B_MJN, to="FIPS") ~ 1,
    fipsst %in% cdlTools::fips(CONTROL_B_MJN, to="FIPS") ~ 0),
  event_tx_balance_fgh = case_when(
    fipsst %in% cdlTools::fips(TREATED_B_FGH, to="FIPS") ~ 1,
    fipsst %in% cdlTools::fips(CONTROL_B_FGH, to="FIPS") ~ 0),
)

# Sad or hopeless
event_sad_1b <- felm(sad_hopeless ~ event_tx_balance_sad*year_cat |
                       year + fipsst + birth_year | 0 | fipsst,
                     data    = yrbs_event,
                     weights = yrbs_event$weight_2)
event_sad_2b <- felm(sad_hopeless ~ event_tx_balance_sad*year_cat +
                       age_2 + sex_2 + race_eth_2 + grade_2 +
                       elig_1_5 + elig_6_18 + has_eitc + federal_pct + refundable + max_bft_3 |
                       year + fipsst + birth_year | 0 | fipsst,
                     data    = yrbs_event,
                     weights = yrbs_event$weight_2)

# Considered suicide
event_con_1b <- felm(cons_suicide ~ event_tx_balance_con*year_cat |
                       year + fipsst + birth_year | 0 | fipsst,
                     data    = yrbs_event,
                     weights = yrbs_event$weight_2)
event_con_2b <- felm(cons_suicide ~ event_tx_balance_con*year_cat +
                       age_2 + sex_2 + race_eth_2 + grade_2 +
                       elig_1_5 + elig_6_18 + has_eitc + federal_pct + refundable + max_bft_3 |
                       year + fipsst + birth_year | 0 | fipsst,
                     data    = yrbs_event,
                     weights = yrbs_event$weight_2)

# Attempted suicide
event_att_1b <- felm(suicide_att ~ event_tx_balance_att*year_cat |
                       year + fipsst + birth_year | 0 | fipsst,
                     data    = yrbs_event,
                     weights = yrbs_event$weight_2)
event_att_2b <- felm(suicide_att ~ event_tx_balance_att*year_cat +
                       age_2 + sex_2 + race_eth_2 + grade_2 +
                       elig_1_5 + elig_6_18 + has_eitc + federal_pct + refundable + max_bft_3 |
                       year + fipsst + birth_year | 0 | fipsst,
                     data    = yrbs_event,
                     weights = yrbs_event$weight_2)

# Alcohol use
event_alc_1b <- felm(alcohol ~ event_tx_balance_alc*year_cat |
                       year + fipsst + birth_year | 0 | fipsst,
                     data    = yrbs_event,
                     weights = yrbs_event$weight_2)
event_alc_2b <- felm(alcohol ~ event_tx_balance_alc*year_cat +
                       age_2 + sex_2 + race_eth_2 + grade_2 +
                       elig_1_5 + elig_6_18 + has_eitc + federal_pct + refundable + max_bft_3 |
                       year + fipsst + birth_year | 0 | fipsst,
                     data    = yrbs_event,
                     weights = yrbs_event$weight_2)

# Marijuana use
event_mjn_1b <- felm(marijuana ~ event_tx_balance_mjn*year_cat |
                       year + fipsst + birth_year | 0 | fipsst,
                     data    = yrbs_event,
                     weights = yrbs_event$weight_2)
event_mjn_2b <- felm(marijuana ~ event_tx_balance_mjn*year_cat +
                       age_2 + sex_2 + race_eth_2 + grade_2 +
                       elig_1_5 + elig_6_18 + has_eitc + federal_pct + refundable + max_bft_3 |
                       year + fipsst + birth_year | 0 | fipsst,
                     data    = yrbs_event,
                     weights = yrbs_event$weight_2)

# Physical fight
event_fgh_1b <- felm(fight ~ event_tx_balance_fgh*year_cat |
                       year + fipsst + birth_year | 0 | fipsst,
                     data    = yrbs_event,
                     weights = yrbs_event$weight_2)
event_fgh_2b <- felm(fight ~ event_tx_balance_fgh*year_cat +
                       age_2 + sex_2 + race_eth_2 + grade_2 +
                       elig_1_5 + elig_6_18 + has_eitc + federal_pct + refundable + max_bft_3 |
                       year + fipsst + birth_year | 0 | fipsst,
                     data    = yrbs_event,
                     weights = yrbs_event$weight_2)

# Get values from models
balance_df <- NULL

# Adolescents (FE only)
balance_df <- make_event_df_2(balance_df, event_sad_1b, "Adolescents (FE only)")
balance_df <- make_event_df_2(balance_df, event_con_1b, "Adolescents (FE only)")
balance_df <- make_event_df_2(balance_df, event_att_1b, "Adolescents (FE only)")
balance_df <- make_event_df_2(balance_df, event_alc_1b, "Adolescents (FE only)")
balance_df <- make_event_df_2(balance_df, event_mjn_1b, "Adolescents (FE only)")
balance_df <- make_event_df_2(balance_df, event_fgh_1b, "Adolescents (FE only)")

# Adolescents (fully adjusted)
balance_df <- make_event_df_2(balance_df, event_sad_2b,
                              "Adolescents (fully adjusted)")
balance_df <- make_event_df_2(balance_df, event_con_2b,
                              "Adolescents (fully adjusted)")
balance_df <- make_event_df_2(balance_df, event_att_2b,
                              "Adolescents (fully adjusted)")
balance_df <- make_event_df_2(balance_df, event_alc_2b,
                              "Adolescents (fully adjusted)")
balance_df <- make_event_df_2(balance_df, event_mjn_2b,
                              "Adolescents (fully adjusted)")
balance_df <- make_event_df_2(balance_df, event_fgh_2b,
                              "Adolescents (fully adjusted)")

# Reorder outcomes
balance_df$outcome <- factor(
  balance_df$outcome, levels=c("Sad or hopeless in past year", "Considered suicide in past year", "Attempted suicide in past year", "Alcohol use in past month", "Marijuana use in past month", "Physical fight in past year"))

# Get min. and max. N
balance_n <-
  (c(length(event_sad_1b$residuals), length(event_sad_2b$residuals),
     length(event_con_1b$residuals), length(event_con_2b$residuals),
     length(event_att_1b$residuals), length(event_att_2b$residuals),
     length(event_alc_1b$residuals), length(event_alc_2b$residuals),
     length(event_mjn_1b$residuals), length(event_mjn_2b$residuals),
     length(event_fgh_1b$residuals), length(event_fgh_2b$residuals)))
min(balance_n); max(balance_n)

# Generate event study plot: Balanced panel
event_plot_balance <- print_event_plot(
  balance_df,
  Y_TITLE    = "Effect of raising minimum wage\non indicated outcome by year",
  Y_MIN      = -0.1,
  Y_MAX      =  0.1
)

# Export figure
ggsave(plot=event_plot_balance, file="Exhibits/YRBS event studies, balanced.pdf",
       width=6, height=6, units='in', dpi=600)
